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INTRODUCTION 


This report is the fourth of a series of reports which summarizes 
the research activities of the Geodynamics Branch. The branch is 
located within the Laboratory for Terrestrial Physics of the Space 
and Earth Sciences Directorate of the Goddard Space Flight Center. 
The branch's major research activites focus on solid-earth dynamics 
including crustal movements, ocean topography and dynamics, and 
planetology. For the solid-earth, research is conducted in 
tectonic plate motion, lithospheric structure and dynamics, regional 
scale crustal deformation, polar motion, earth orientation and 
rotation, and gravity field modeling. Satellite laser ranging and 
very long baseline interferometric data provide much of the obser- 
vational information used in studying these topics. For the oceans, 
research is conducted in ocean topography and circulation, and 
mesoscale variability and structure. Analysis of satellite alti- 
meter data, ocean geoid modeling, and precision orbit determination 
are fundamental for performing this research. Planetology studies 
are of recent interest to the branch and consist of geodetic and 
geodynamic analyses. Impetus was given to this research by the 
future Mars Observer Space Mission in which the branch plans to 
play a key role. 

The NASA programs which are supported by the work described in this 
document, include the Geodynamics and Ocean Programs, the Crustal 
Dynamics Project, the Ocean Topographic Experiment (TOPEX) and the 
Mars Observer Mission. The report highlights the investigations 
conducted by the Geodynamics Branch Staff and its on-site contrac- 
tors during calendar year 1985. The individual papers, are grouped 
into the chapters: Tectonophysics, Sea Surface Topography, Geodesy 
and Geopotential Field, Planetology, and Geophysical Techniques. 
More detailed information on the activities of the branch or the 
particular research efforts described herein, can be obtained 
through the branch office or from individual staff members. 
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CHAPTER I 


TECTONOPHYSICS 


OVERVIEW 


For the past two decades, the theory of plate tectonics has pro- 
vided a framework for studying a variety of geodynamic processes. 
Foremost among these is the motion of the plates themselves. 
Measurements using space techniques are particularly well -suited to 
the direct determination of tectonic plate kinematics. In fact, 
measurements using satellite laser ranging and very long baseline 
interferometry have matured to the point where three-dimensional 
descriptions of individual site locations are being determined. 
This new information supplements the interesting, but more limited, 
determinations of intersite baselines which have been obtained for 
several years. Currently, changes in site positions and baseline 
lengths are being studied from sites situated on most of the major 
tectonic plates. This year's reports concentrate on measurements 
involving plate motion and deformation in the western United States 
and the Pacific. Particularly interesting are the comparisons 
between the currently observed rates and those deduced from the 
geological record. In many cases the contemporary and geologic 
rates closely agree, but the remaining differences are quite 
intriguing. The continued monitoring of these motions will address 
the questions of whether temporal variations exist in the motions, 
how deformations occur near plate boundaries, and whether current 
models of plate motion are adequate. The kinematic descrip- 
tions of the tectonic plate motion provides input for dynamical 
models of the associated driving processes. Howevever, the kine- 
matic description of plate tectonics can be determined directly by 
space measurements, while the determination of driving forces asso- 
ciated with geodynamical processes requires a merging of observa- 
tional data and geophysical models. Accordingly several studies 
conducted within the Geodynamics Branch involve analytical and 
numerical modeling of the dynamics of tectonic phenomena. During 
the past year, models have been developed of the lithospheric 
deformation associated with continental collisions, the deformation 
associated with folding in the Indian Ocean, the stresses asso- 
ciated with the Tibetan and Himalayan uplifts, and the changes in 
earth rotation, polar motion and the gravity field induced by 
earthquakes. 
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The dynamical processes that are responsible for some of the 
tectonophysics phenomena discussed in this section are intimately 
related to those discussed in subsequent chapters of this report on 
Geodesy and Geopotential and on Planetology. The interested reader 
is referred to these sections for additional reports on the earth 
and planetary physics studies being conducted within the Geodynamics 
Branch. 

Contributors to this chapter are: B. Fong Chao, Demosthenes C. 
Chri stodoul i di s , Thomas A. Clark, Steven C. Cohen, David Gordon, 
Richard S. Gross, Han-Shou Liu, Chopo Ma, James W. Ryan, David E. 
Smith, and Maria T. Zuber. 
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TECTONIC MOTIONS IN WESTERN USA FROM SLR 


Demosthenes C. Christodoul idis 
David E. Smith 


OBJECTIVE 

The objective of this work is to derive tectonic motion vectors for 
tracking sites in the western United States from analysis of LAGEOS 
satellite laser tracking data. 


BACKGROUND 

The analysis of laser satellite tracking data utilizes principles 
from dynamic satellite geodesy. Products of this type of analysis 
are precise three-dimensional laser tracking coordinates. 
Dynamic analyses have been performed on ,data accumulated over a 6- 
year period (1979-1984) resulting in precise annual tracking station 
positions. During this 6-year period, laser systems and data 
consistency have reached advanced levels of maturity yielding a 
strong network geometry and better network resolution than in 
previous years. The analysis of the data at GSFC yielded the SL6 
solution consisting of improved Earth rotation, station coordinates 
and geodetic constants. In this work, the original SL6 solution 
was iterated with the introduction of better nutation and ocean 
tide models. Annual geodesic lengths were computed from the precise 
final annual station positions in our new SL6 system. 


RECENT ACCOMPLISHMENTS 

Absolute velocity vectors for the laser tracking stations in the 
western United States have been computed via a least-squares 
adjustment model. "Absolute" in this context means that the veloci- 
ties are referenced with respect to mantle intrusions into the 
Earth's crust, known as hotspots. In this model, two types of 

stations are defined: "external" stations which have their abso- 

lute plate motion defined by a geologic plate motion model (these 
are shown in Figure 1); and "internal" stations for which the 
absolute motion is sought (shown in Figure 2). The Minster and 
Jordan (1978) AMI-2 model was adopted to describe the absolute 
motion of the external stations for two reasons. First, it agrees 
well with the SLR observed rates for these stations, and secondly 
it provides the means to relate the computed velocity vectors to 
the mantle-based reference frame. All the "external" stations 
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have dense tracking data sets and the SLR derived relative rates for 
these stations are within one standard deviation from the relative 
rates implied by the geologic model. The methodology developed in 
this work begins by fitting a slope to the time series of geodesics 
between all pairs of stations which have in common at least two 
annual geodesic lengths. These slopes, or geodesic rates, assumed 
constant are then processed simultaneously in a least-squares 
algorithm to yield absolute velocity vectors for the "internal" 
station network. The errors in the station positions, as approxi- 
mated by the error in the height of the station, are propagated 
throughout the system. The uncertainties in the geodesic rates are 
used as "observational" weights in the simultaneous solution. 

Details of this work are described in Christodoulidis et al . (1986). 


SIGNIFICANCE 

Satellite laser ranging has advanced to the point of providing the 
capability to monitor tectonic motions. The analysis made in this 
work guarantees self-consistency in the solution for absolute motion 
since all of the constructable geodesic rates enter into the set of 
simultaneous equations. This is a considerable improvement over 
the previous methods of analyzing the results in a baseline-by- 
baseline manner. The resulting tectonic motion model, the SL6 
Tectonic Model, has been compared with the geologic models of 
Minster and Jordan (1978) and Chase (1978). The SL6 computed 
absolute motions are listed in Table 1 along with motions implied 
by the geologic models. The modeled geodesic rates which traverse 
the San Andreas Fault indicate that the rate of change (shortening) 
over these lines has slowed considerably in the last 3 years. The 
best estimate for the period 1982-84 is -2.56 cm/year. The slower 
northwesterly absolute motion for the Monument Peak site is largely 
responsible for this behavior. The current results indicate that 
the motion along the SAFE line (San Andreas Fault Experiment line, 
Momument Peak [7110] to Quincy [7051 & 7109] now agrees with the 
rate obtained from the local ground survey networks. This signifi- 
cant change in the rate of the SAFE line raises the question 
whether Monument Peak is located on the Pacific Plate or is on some 
California platelet. Implementation of dense regional observation 
networks and continuous analysis of SLR and VLBI data are necessary 
for furthering our understanding of the dynamic behavior of this 
enti re region. 


FUTURE EMPHASIS 

The accuracy -and reliability of SLR Tectonic Motion Models will 
be improved as additional data are acquired and processed. The 
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least-squares algorithm is being improved to propagate full covar- 
iance information and the analysis portions are being supplemented 
to better locate aberrant data. With these improvements, it is 
expected that these results will have considerable impact among the 
tectonophysics community. 


REFERENCES 

Chase, C.6., "Plate Kinematics: The Americas, East Africa, and the 

Rest of the World," Earth & Plan. Sci . Letters, Vol . 37. dd 355- 
368,1978. 

Christodoulidis, D.C., D.E. Smith, S.M. Klosko and J.W. Robbins, 
"Tectonic Motion in Western U.S.A. from Satellite Laser Ranging," 
Proc. of AAPG Memoir Workshop, Corvallis, Oregon, in press, 1986. 

Minster, J.B. and T.H. Jordan, "Present-Day Plate Motions," 
J. Geophys. Res. . Vol. 83, pp 5331-5354, 1978. 
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Table 1 . A Comparison of Absoluted Station Motions From Differing Models 


7210 

7090 

7907 

7105 

7086 

7051 

7062 

7109 

7110 
7112 
7114 
7122 


MAGNITUDE (CM/YR) 


MINJOR 

CHASE 

SL6 

9.68 

9.10 

9.68 

7.97 

7.24 

7.97 

3.12 

2.34 

3.12 

2.69 

2.71 

2.69 

2.67 

2.69 

2.95 

2.41 

2.79 

3.86 

6.38 

6.50 

4.98 

2.41 

2.79 

1.19 

6.31 

6.44 

4.50 

2.55 

2.77 

3.12 

2.49 

2.79 

3.30 

2.72 

2.61 

4.15 


AZIMUTH (DEG) 


MINJOR 

CHASE 

SL6 

300.40 

297.14 

300.40 

21.71 

29.11 

21.71 

265.76 

288.42 

265.76 

251.53 

264.95 

251.53 

241.14 

242.82 

239.86 

233.86 

232.01 

229.63 

296.64 

296.37 

307.53 

233.86 

232.01 

251.66 

296.62 

296.38 

301.14 

239.39 

243.17 

346.89 

235.46 

233.66 

174.31 

241.06 

239.91 

201.74 
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Figure 2. “Internal” SLR network. 
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DETERMINATION OF SITE MOTIONS IN THE WESTERN UNITED STATES FROM 

MARK III VLBI MEASUREMENTS 


David Gordon 
Thomas A. Clark 


OBJECTIVE 

The objective of the work reported here is to determine site 
motions at 14 fixed and mobile VLBI sites in the western continental 
United States using a large data set of Mark III VLBI measurements. 


BACKGROUND 

For several years NASA's Crustal Dynamics Project (CDP) has used 
the Mark III VLBI system to study crustal deformations and movements 
in the western United States. Several dozen sites, mostly in 
California, have been visited repeatedly by 3 mobile VLBI units. 
Additionally, 4 fixed sites in this region have been repeatedly used 
in conjunction with the mobile sites. The CDP is currently per- 
forming some 2 dozen experiments yearly involving these sites. 
VLBI data involving some 3 years of mobile and some 5 years of 
fixed station baseline measurements to sites in the western U.S. 
have been analyzed. Sufficient data now exist to begin assembling 
a detailed picture of the plate motions in this region. 

For a VLBI geodetic experiment, the quantity of interest is the 
baseline vector between the two antennas and how this vector changes 
with time. (For the mobile systems, the baseline vector is trans- 
posed to a fixed geodetic marker). The standard technique is to 
define an a priori baseline, as measured at some reference epoch, 
and to monitor changes in that baseline each time it is remeasured. 
Baseline changes can be broken up into three orthogonal components: 
length, transverse and vertical. The rates of change in these 
three directions are simply the length, transverse, and vertical 
velocities. Length changes are simply increases (positive) or 
decreases (negative) in the length of the baseline. The transverse 
component is perpendicular to the baseline vector and directed 
along the horizon at either site. Positive transverse motion is 
defined as a clockwise rotation of the baseline as seen from above. 
The vertical component is perpendicular to both the length and 
transverse components, and is directed towards the vertical at the 
midpoint of the baseline. A positive vertical change means either 
an increase in height at station §1 or a decrease at station #2. 

preceding page ?lank not HU.ED 
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Measurements of transverse and vertical components are strongly 
affected by errors in the determination of the Earth's pole position 
and UT1. The vertical component also suffers strongly by uncertain- 
ties in the tropospheric delay calibration. The length component, 
however, is essentially unaffected by errors in the Earth orienta- 
tion, and length errors due to tropospheric delay calibration 
errors can generally be held small (to about the 1 cm level for the 
baseline lengths studied here) through the consistent use of sur- 
face meteorology measurements. Until recently, only the length 
components of baseline velocities have been useful in attempts at 
resolving individual site motions. The use of the transverse 
component requires the availability of a highly accurate and homo- 
geneous Earth orientation series. Until recently, no such series 
were available; BIH Circular D is not accurate enough and other, 
more accurate series were not homogeneous. 


RECENT ACCOMPLISHMENTS 

Since late December, 1983, the National Geodetic Survey has been 
conducting program IRIS, utilizing a 4 station (Westford, Ft. Davis, 
Richmond and Wettzelll) VLBI network to very accurately determine 
pole position and ujri at 5-day intervals. The GSFC VLBI group has 
compiled and regularly updates an Earth orientation series based on 
these IRIS measurements. (NGS also maintains an IRIS Earth 
orientation series. The two series differ slightly mainly in the 
use of different a priori i site and source positions). This series, 
which is continuous since the beginning of 1984, is homogeneous and 
appears to be highly accurate. During the past year, software has 
been developed to map the Earth orientation parameters from this 
IRIS series into data bases at the time of analysis. Using this 
technique, errors in the transverse baseline component are limited 
to the 1 - 2 cm level for baselines in the western U.S. This has 
allowed accurate determination of transverse velocities using data 
from 1984 and 1985. 

A data set comprised of baseline measurements between 14 mobile and 
fixed VLBI sites in the western continental U.S. plus measurements 
from Westford and Gilmore Creek (Fairbanks) to these sites has been 
compiled for this study. This data set contains some 500 individual 
baseline measurements (not counting several hundred measurements of 
Westford to Ft. Davis as part of the IRIS and POLARIS programs) on 
some 57 individual baselines. For each baseline, a length and 
transverse velocity was computed using a simple linear least squares 
fitting program. The transverse velocities were computed using 
only data from 1984 and 1985 (about 360 individual baseline 
measurements) since the IRIS series does not exist before this. 
Table 1 lists, for each baseline, the length and transverse veloci- 
ties in mm/year along with their computed standard deviation. In 
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order to get a more realistic estimate of the standard deviation, 1 
cm of noise was added to each length measurement and 2 cm of noise 
was added to each transverse measurement. As an example of the 
data, in Figure 1 plots of the length and transverse data for the 
Mojave to Monument Peak baseline are shown. 

Software has been developed to perform a least squares analysis of 
this set of length and transverse velocities and determine indivi- 
dual site motions. In the analysis performed, the velocities were 
weighted by the inverse of their standard deviations. It was 
necessary to hold two stations fixed to define a reference frame. 
Westford, included because it helps to resolve the motion at Ft. 
Davis, was held fixed as one of the anchor points. OVRO was chosen 
as the other anchor point; its motion was restricted to less than 1 
mm/year. OVRO was chosen because it is well east of the San Andreas 
fault and because there exists a long history of VLBI between 
Westford and OVRO showing no significant motion. Fairbanks was 
included because it helps to resolve the northward motions in 
California. The program allowed specifying a priori first guess 
site motions as well as restricting a site to a specified a priori 
value with varying degrees of rigidity. The final solution obtained 
is thus not unique. However, various solutions with different 
constraints seem to differ at only the few mm/year level in terms 
of site motions. Table 2 gives the results of the final solution. 
Site speed in mm/year and azimuthal direction in degrees are given. 
Graphically, the results are presented on a map of the western U.S. 
in Figure 2 (Clark, 1986). The arrows from each site indicate its 
direction of motion and relative speed. Figure 2 also shows error 
ellipses. These are estimates based on an examination of the data 
and are not computed by the least squares fitting program. 


SIGNIFICANCE 

This work yields the most detailed picture yet of crustal motions 
in the western U.S. from space geodetic techniques. Seven sites 
situated near or west of the San Andreas fault show large, signifi- 
cant motions. They are Monument Peak, Vandenberg, JPL (Pasadena), 
Pearblossom, Fort Ord, Presidio (San Francisco), and Pt. Reyes. 
While a detailed analysis of these motions is not attempted here, a 
few points are worth noting. The direction of motion as one moves 
northward through these sites (except Presidio) rotates to follow 
the direction of the San Andreas fault. The Presidio site, in San 
Francisco, appears anomalous. It is just east of the San Andreas 
and shows westward motion towards the fault, motion which is not 
understood at this time. At Monument Peak, only about 60% of the 
full motion of the Pacific Plate relative to the North American is 
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seen, indicating that the other 40% is taken up in offshore faults. 
However, Vandenberg and Fort Ord show motion of some 4-6 cm/year, 
which accounts for most of, if not all, the Pacific Plate motion 
and indicates that little if any slippage is occurring in the faults 
offshore of these regions. Thus some 4-6 cm/year of slippage 
and/or strain accumulation must be occurring in the region from 
Vandenberg/Fort Ord eastward to the San Andreas. This could be 
occurring among the many faults of the Hosgri and San Gregorio 
fault systems or along the San Andreas or a combination of these. 

Additionally, Ft. Davis shows clear eastward motion of about 7 
mm/year, in general agreement with earlier findings. There seems 
to be motion at Platteville and Quincy but it is very uncertain at 
this time. Also, the southward motion of Fairbanks relative to 
eastern California is not considered significant. 


FUTURE EMPHASIS 

During the next year, additional measurements are planned to 
Vandenberg, Monument Peak, Yuma, Platteville, JPL, and Quincy, as 
well as the fixed sites Hatcreek, OVRQ, Mojave, and Ft. Davis. 
These can be expected to strengthen the site motion determinations 
at these sites. Also four sites which have been visited before, 
Pinyon Flats, Flagstaff, Ely and Blackbutte will be revisited and 
enough data should then be obtained to include them in a site motion 
solution. Additional measurements from Fairbanks to Mojave, 0VR0 
and Hatcreek in 1986 should help to clarify the relation between 
Alaska and eastern California. 


REFERENCE 

Clark, Thomas A., "Plate Deformation Measurements Using VLBI in the 
Western U.S.," Report presented to the CDP Investigators Meeting, 
JPL, Pasadena, California, March 27, 1986. 
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Table 1. Baseline Motion Information 


Baseline 


Rate of Change of 
Baseline Length 

Number of 
Data Sets 

Solution 

Residual 

Normalized 

Solution 

Residual 

Mojave 

to Vandenberg 

(mm/year) 
Value Sigma 

37 12 

8 

(mm/year) 

17 

(unitless) 

1.4 

Mojave 

to Gilcreek 

-20 

10 

17 

-2 

0.2 

Mojave 

to Kauai 

28 

11 

13 

6 

0.5 

Mojave 

to Kwajalein 

-14 

28 

11 

-6 

0.2 

Mojave 

to Kashima 

-5 

25 

15 

-5 

0.2 

Vandenberg 

to Gilcreek 

-72 

13 

16 

17 

1.3 

Vandenberg to Kauai 

-7 

16 

9 

-5 

0.3 

Gilcreek 

to Kauai 

-34 

16 

12 

3 

0.2 

Gilcreek 

to Kwajalein 

-35 

16 

10 

8 

0.5 

Gilcreek 

to Kashima 

-7 

20 

12 

-19 

1.0 

Kauai 

to Kwajalein 

-2 

17 

11 

30 

1.8 

Kauai 

to Kashima 

-69 

25 

10 

-15 

0.6 

Kwajalein 

to Kashima 

-87 

17 

8 

-17 

1.0 


Table 2. Site Motion Solution 



Motion 

Azimuth 

STATION 

(mm/yr.) 

(degrees) 

MOJAVE 12 

0 . 

— 

KASHIMA 

0 . 

— 

VNDNBERG 

62. 

326. 

GILCREEK 

19. 

136. 

KAUAI 

61. 

308. 

KWAJAL26 

76. 

333. 


! 
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TRANSVERSE RESIDUAL (cm) LENGTH - 27405584.392 cm 


MOJAVE TO MONUMENT PEAK BASELINE LENGTH MEASUREMENTS 




Figure 1. Baseline length and transverse data versus time for the baseline from Mojave to Monument 
Peak as measured by Mark III VLBI. The vertical line for each point represents the error 
bar for each individual baseline measurement, where 1 cm of noise was added to each 
length standard deviation and 2 cm of noise was added to each transverse standard 
deviation. The least squares fit and its slope for each plot represent the length and 
transverse velocities of this baseline. 
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Figure 2. Relative site motions and their estimated error ellipses for VLBI sites in the western U.S. 






PACIFIC PLATE MOTION OBSERVED WITH MARK III VLBI 


James W. Ryan 

! Chopo Ma 

i 


OBJECTIVE 

A primary objective of the NASA Crustal Dynamics Project (CDP) is 
to make contemporary measurements of the large-scale motions of the 
Earth's crust. These data are then the basis from which models of 
the current dynamics of the Earth can be made. Because of its 
seismicity, one area of great interest is the rim of the Pacific. 


BACKGROUND 

As a part of its program to monitor global -scale tectonic plate 
motions the CDP in early 1984 began a series of VLBI experiments to 
measure the motions of sites on the Pacific plate relative to sites 
on the Pacific rim. The Pacific plate sites, called here Kauai, 
Kwajalein, and Vandenberg, are located in the Hawaiian and Marshall 
Islands (3700 kilometers west Hawaii) and at Point Arguello in 
California, respectively. The Pacific rim sites, called Mojave, 

Gil creek, and Kashima, are located in California's Mojave desert, 
in central Alaska, and in Japan, respectively. 

The first two experiments in the Pacific motion series took place 
in January and February 1984 between Mojave and Kashima and were 
largely engineering tests of the new Japanese-built K-3 VLBI data 
acquisition system. In the summers of 1984 and 1985 there were two 
series of six experiments each using various subsets of the sites 
mentioned above. Vandenberg also participated with Gil creek in a 
series of Alaskan regional experiments in 1984 and 1985. Table 1 
shows the number of data sets for the various baselines. Some 
experiments produced two data sets. 


RECENT ACCOMPLISHMENTS 

These data were analyzed as a part of a large least-squares adjust- 
ment of all the CDP and IRIS Mark III VLBI data acquired in 1984 
and 1985. (IRIS is the operational Earth rotation monitoring 
program of the National Geodetic Survey which uses Mark III VLBI). 
In all there were 167 data sets. The parameters estimated in this 
solution were experiment-dependent site positions, globally 
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estimated radio source coordinates, experiment-dependent offsets to 
the nutations in obliquity and longitude, and time and station- 
dependent clock and atmosphere parameters. Since the site coordi- 
nates were estimated in this solution for each epoch independently, 
any tectonic effects which might have occurred during the 2-year 
period could be accommodated. The results of this solution were 
condensed in the form of tables of baseline length during the time 
interval. Errors in the a priori values of pole position and UT1 
from the B.I.H. do not affect these length histories since they 
affect only the orientation of the networks in the VLBI coordinate 
system. From each of these baseline histories a length and rate of 
change of length were determined. Table 1 shows the rates of 
change of baseline length with their formal standard error from the 
least-squares fits. The rates of change reflect possible plate 
motion effects. 

In order to test whether these observed motions are consistent with 
a simple model of Pacific plate motion relative to the North American 
plate, the 13 rates-of-change given in Table 1 were used to estimate 
by weighted least-squares the motions of Kwajalein, Kauai, 
Vandenberg, and Gilcreek. For each station the local north and east 
components of the station velocity were estimated in a coordinate 
frame in which Mojave and Kashima were held fixed. This is an ad hoc 
constraint. However, there are valid geophysical reasons (not 
discussed here) for assuming that Kashima and Mojave have no rela- 
tive motion, and the direct baseline measurements between these 
sites are consistent with no motion. In any case some such con- 
straint is needed to anchor the coordinate frame. Table 2 presents 
the results of this estimation as the total station motion in the 
horizontal plate and the azimuth of the motion. The uncertainties 
are not tabulated. Since this data set spans only 1.5 years with 
only 13 data points to estimate 6 parameters, the uncertainties 
produced by the estimator are not meaningful. The justification for 
the correctness of this solution is intuitive and rests on its 
ability to fit the observed motions. 


SIGNIFICANCE 

The last two columns in Table 1 present the differences between the 
observed and computed motions and these differences normalized by 
the corresponding rate-of-change sigma. Three baseline velocities, 
Mojave to Vandenberg, Vandenberg to Gilcreek and Kauai to Kwajalein, 
have large normalized residuals indicating that these velocities 
and sigmas are not entirely consistent with the other data. Except 
for these three values the residuals are 1 sigma or less and eight 
are at the level of 0.5 sigma or less. Table 2 shows that 
Vandenberg, Kauai and Kwajalein are all moving with a velocity of 
60-70 mm/yeaF in a northwesterly direction. This is entirely 
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consistent with current models for the motion between the Pacific 
and North American plates. The solution also indicates a relative 
motion of approximately 15 mm/year between Kauai and Kwajalein and 
motion of Gilcreek towards Mojave at a rate of 19 mm/year. Neither 
of these smaller motions is consistent with current plate motion 
models. 

In summary, a motion of the Pacific plate relative to the North 
American plate of approximately 65 mm/year in a northwesterly 
direction appears to have been detected. The solution also indi- 
cates that there are intraplate motions at the level of 20 mm/year, 
which are at variance with current plate motion models. 


FUTURE EMPHASIS 

Experiments involving the current Pacific stations planned for the 
summer of 1986 will nearly double the time span of the data sets 
on these baselines. This will significantly improve the precision 
of the observed motions and provide the basis for a statistically 
significant estimation of the motion of the stations. During a 
trip made by the second author to China and Japan in 1985, discus- 
sions were held with researchers at the Shanghai Observatory, the 
Yunnan Observatory, the Research Institute of Mappinq and Surveying 
in Beijing, and the Kashima Space Research Center about future VLBI 
programs which could extend the baselines around the Pacific into 
the Asian mainland. An exploratory VLBI experiment has been per- 
formed between Shanghai and Kashima. The Shanghai Observatory has 
acquired Mark III VLBI equipment, and its new 25-m antenna is 
planned to be operational by 1986. 


Table 1. 


Baseline 

WESTFORD — FT. DAVIS 
WESTFORD - OVRO 
WESTFORD — FAIRBANKS 
FT. DAVIS - PLATTVIL 
FT. DAVIS — YUMA 
FT. DAVIS — MON PEAK 
FT. DAVIS - MOJAVE 
FT. DAVIS - OVRO 
FT. DAVIS — VNDNBERG 
FT. DAVIS - QUINCY 
FT. DAVIS - HATCREEK 
FT. DAVIS - FAIRBANKS 
PLATTVIL - MOJAVE 
PLATTVIL - OVRO 
PLATTVIL - HATCREEK 
YUMA - MON PEAK 
YUMA — MOJAVE 
YUMA - OVRO 
YUMA - VNDNBERG 
MON PEAK — MOJAVE 
MON PEAK - OVRO 
MON PEAK - VNDNBERG 
MON PEAK — QUINCY 
MON PEAK - HATCREEK 
MOJAVE — PBLOSSOM 
MOJAVE - PASADENA 
MOJAVE — OVRO 
MOJAVE — VNDNBERG 
MOJAVE - QUINCY 
MOJAVE - HATCREEK 
MOJAVE - FORT ORD 
MOJAVE - PRESIDIO 
MOJAVE - PT REYES 
MOJAVE - FAIRBANKS 
PBLOSSOM - PASADENA 
PBLOSSOM - OVRO 
PBLOSSOM - VNDNBERG 
PASADENA - OVRO 
PASADENA - VNDNBERG 
OVRO - VNDNBERG 
OVRO — QUINCY 
OVRO - HATCREEK 
OVRO — FORT ORD 
OVRO — PRESIDIO 
OVRO — PT REYES 
VNDNBERG - QUINCY 


LENGTH 

VELOCITY 

mm/yr 


- 9 . 

+ /- 

2 . 

2 . 

+/- 

2 . 

5 . 

+ /- 

5 . 

- 13 . 

+/- 

8 . 

11 . 

+/- 

4 . 

39 . 

+ /- 

5 . 

5 . 

+/- 

2 . 

8 . 

+ /- 

2 . 

47 . 

+/- 

7 . 

- 0 . 

+/- 

9 . 

15 . 

+/- 

4 . 

0 . 

+/- 

22 . 

- 9 . 

+/- 

8 . 

- 1 . 

+/- 

8 . 

- 3 . 

+/- 

7 . 

32 . 

+/- 

8 . 

5 . 

+ /- 

5 . 

8 . 

+/- 

10 . 

44 . 

+/- 

11 . 

- 25 . 

+/- 

4 . 

- 31 . 

+/- 

4 . 

15 . 

+ /- 

6 . 

- 35 . 

+/- 

7 . 

- 22 . 

+ /- 

6 . 

9 . 

+/- 

3 . 

10 . 

+ /- 

6 . 

- 2 . 

+/- 

2 . 

18 . 

+ /- 

3 . 

- 2 . 

+/- 

4 . 

5 . 

+/- 

3 . 

31 . 

+/- 

6 . 

28 . 

+/- 

4 . 

29 . 

+ /- 

6 . 

- 10 . 

+/- 

5 . 

- 2 . 

+/- 

6 . 

- 7 . 

+ /- 

3 . 

13 . 

+/- 

4 . 

- 17 . 

+ /- 

6 . 

19 . 

+/- 

10 . 

- 10 . 

+/- 

3 . 

- 4 . 

+/- 

4 . 

5 . 

+/- 

3 . 

13 . 

+ /- 

4 . 

25 . 

+/- 

6 . 

29 . 

+/- 

- 6 . 

- 44 . 

+ /- 

- 17 . 


TRANSVERSE 

VELOCITY 

mm/yr 


- 7 . +/- 10 . 
0 . +/- 100 
20 . +/- 10 . 
7 . +/- 4 . 
4 . +/- 5 . 
55 . +/- 16 . 
37 . +/- 10 . 
2 . +/- 10 . 
31 . +/- 20 . 
10 . +/- 10 . 
23 . +/- 10 . 
.6 +/- 10 . 
16 . +/- 100 
1 . +/- 8 . 
0 . +/- 100 

19 . +/- 13 . 
29 . +/- 5 . 
22 . +/- 5 . 

34 . +/- 7 . 

36 . +/- 22 . 

20 . +/- 6 . 

27 . +/- 5 . 

31 . +/- 6 . 
- 0 . +/- 2 . 
29 . +/- 4 . 
19 . +/- 6 . 

2 . +/- 4 . 

35 . +/- 6 . 
- 11 . +/- 5 . 

22 . +/- 7 . 
0 . +/- 10 . 
- 4 . +/- 10 . 
19 . +/- 5 . 
7 . +/- 8 . 

19 . +/- 8 . 

- 6 . +/- 10 . 
53 . +/- 4 . 
19 . +/- 7 . 
5 . +/- 4 . 
53 . +/- 8 . 
10 . +/- 7 . 

39 . +/- 8 . 
52 . +/- 28 . 
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Table 1. (Con’t) 




LENGTH 

TRANSVERSE 



VELOCITY 

VELOCITY 

Baseline 

mm/yr 

mm/yr 

VNDNBERG 

- HATCREEK 

-28. +/- 4. 

34. +/- 5. 

VNDNBERG 

- FORT ORD 

8. + /- 5. 

22. +/- 6. 

VNDNBERG 

- PRESIDIO 

-1. +/- 6. 

-23. +/- 7. 

VNDNBERG 

- PT REYES 

2. +/- 9. 

9. +/- 8. 

VNDNBERG 

— FAIRBANKS 

-70.. +/- 6. 

10. +/- 20. 

QUINCY 

- HATCREEK 

3. +/- 10. 

-11. +/- 6. 

HATCREEK 

- FORT ORD 

-25. +/- 5. 

34. +/- 6. 

HATCREEK 

- PRESIDIO 

33. +/- 8. 

39. +/- 9. 

HATCREEK 

- PT REYES 

-23. +/- 8. 

18. +/- 9. 

FORT ORD 

- PRESIDIO 

-9. +/- 11. 

-40. +/- 115 

PRESIDIO 

- PT REYES 

-4. +/- 9. 

23. +/- 20. 
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Table 2. 



MOTION 

AZIMUTH 

STATION 

(mm/yr) 

(degrees) 

WESTFORD 

0 . 

0 . 

FT. DAVIS 

7. 

132. 

PLATTVILLE 

13. 

183. 

YUMA 

2. 

71. 

MONUMENT PEAK 

35. 

304. 

MOJAVE 

2. 

359. 

PEARBLOSSOM 

24. 

308. 

JPL (PASADENA) 

27. 

318. 

OVRO 

0 . 

210. 

VANDENBERG 

50. 

325. 

QUINCY 

19. 

61. 

HATCREEK 

5. 

9. 

FORT ORD 

48. 

334. 

PRESIDIO (SFO) 

36. 

288. 

PT. REYES 

40. 

338. 

FAIRBANKS 

8. 

176. 
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CHANGES IN THE EARTH'S ROTATION AND LOW-DEGREE GRAVITATIONAL FIELD 

INDUCED BY EARTHQUAKES 


B. Fong Chao 
Richard S. Gross 


OBJECTIVE 

The two objectives of this work are: (1) To compute the changes in 
the Earth's rotation parameters (polar motion and length-of-day ) 
and in low-degree gravitational field caused by major earthquakes 
during the period 1977-1985, and (2) To conduct statistical 
analyses on these changes to reveal their geophysical significance. 


BACKGROUND 

The Earth's gravitational field is not static. Instead, it is 
constantly subject to minute changes caused by geophysical varia- 
tions in the density distribution of the Earth's material. These 
density variations also cause the Earth's rotation to fluctuate-- 
they excite the polar motion and change the length-of-day (LOD). 

This study investigates the geodetic and gravitational effect of 
earthquakes as a mechanism by which the Earth changes its density 
distribution. The study of gravitational effect of earthquakes has 
been restricted only to the polar motion, although it has been 
studied many times over in the past (e.g. Lambeck, 1980). No 
conclusive evidence of a significant effect has yet been found. 
Two reasons stand out: (1) inadequacies in accuracy and temporal 

resolution of the conventional polar motion data, and (2) lack of 
truly great earthquakes since the 1964 Alaskan event. However, 
modern techniques, such as satellite laser ranging and very-long- 
baseline-interferometry, have demonstrated the ability to yield 
highly accurate polar motion and LOD data, as well as the ability to 
detect changes in J 2 (e.g. Rubincam, 1984) and other gravitational 
coefficients. While these techniques are improving, future space 
geodetic missions (such as the proposed Geopotential Research 
Mission) may prove to be useful tools to observe the temporal 
changes in the gravitational field (Wagner and McAdoo, 1986). 


RECENT ACCOMPLISHMENTS 

The theoretical formulation is achieved using a normal mode summa- 
tion scheme based on Gilbert (1973) and Gilbert and Dziewonski 
(1975). [The same scheme has been employed by O'Connell and 
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Oziewonski (1976) in computing seismi cal ly-i nduced changes in polar 
motion and LOD, but no complete account of their procedure has been 
published.] Thus we derived our basic formula which turns out to be 
particularly simple both in a mathematical and computational sense: 


ACji,m + = E I n £ M : e n ^ m (r 0 ) 


( 1 ) 


where the symbol : indicates a tensor contraction. 

The left-hand side is the change in the gravitational-field harmonic 
coefficient (the "Stokes Coefficient") of degree i and order m. 
On the right-hand side, I nJl is an integral involving the spheroidal 
ei genel ements of the (spherically symmetric) Earth with overtone 
number n and degree i ; M is the moment tensor of the earthquake; 
and e nJlni (r 0 ) is the elastic strain tensor at the source r 0 asso- 
ciated with the (n, i, m)th spheroidal normal mode. The eigen- 
elements and the strain tensor are obtained from the earth model 
1066B from Gilbert and Oziewonski (1975). The seismic moment 
tensors M are taken from those published regularly by the Harvard 
Group since 1977 (some 2000 events, see Oziewonski et al., 1985). 
The summation in (1) (over n) is found to converge rapidly. We 
have thus computed the seismi cal ly induced changes in some low- 
degree Stokes coefficients. Of particular interest are the 
(m=0) coefficients and the polar motion excitation (t=2, m=l). 

In addition, a formula similar and complementary to (1) gives the 
changes in the length-of-day. 

Major findings are: (1) The earthquake-induced changes we computed 

are in general two orders of magnitude smaller than the observed 
values that are available. (2) Most of these changes show strong 
evidences of non-randomness either in their polarity or in their 
directions. (3) The parameters that show the strongest non-random- 
ness are the dynamic oblateness J2, the total principal moment of 
inertia, the length-of-day, and the sum of as well as the diffe- 
rence between the two equatorial principal moments of inertia. 
They are all inclined toward negative changes, indicating the 
tendency of earthquakes to make the Earth rounder, and to pull in 
mass toward the center of the Earth. (4) The earthquakes continue 
to make the rotational pole drift toward a preferred direction of 
~150°E. This direction is roughly the opposite to that inferred 
from observations. (5) The earthquake-induced changes in the polar 
principal moment of inertia are considerably larger than those in 
its two equatorial counterparts combined. (6) The earthquake- 

induced changes in the zonal harmonics J3, J4, J5 (and perhaps 
higher terms) are much smaller than that in J2; and they show little 
sign of non-randomness. (7) The above behaviors appear to be 
independent of time and the size of the earthquake causing the 
changes. 
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Further investigation will be made on the statistics of the seismic 
excitation of polar motion. The results will be interpreted in the 
hope that light can be shed on the Earth's internal dynamics. On 
the other hand, any great earthquakes in the future will be subjected 
to the study outlined here on an individual basis. The predicted 
geodetic and gravitational changes will be compared to observations 
in search of any possible correlation. 
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FEATURES OF COMPRESSIONAL STRAINING IN CONTINENTAL COLLISIONS 


Steven C. Cohen 


OBJECTIVE 

Continental collisions result in considerable intraplate deforma- 
tion, a geophysical phenomenon that is not wel 1 -incorporated into 
the rigid plate formulation of plate tectonics. The manner in 
which strain and stress are transmitted into the plate interior 
from the plate boundary is an important aspect of this phenomenon. 
The objective of this work is to develop numerical models that 
explain the intraplate deformation and strain propagation in terms 
of the rheological properties of the earth and mechanisms of plate 
i nteracti ons . 


BACKGROUND 

Among the many manifestations of continental collisions are a 
variety of lithospheric and crustal deformations that can extend 
very deep into continental interiors. Mountain building, plateau 
uplift, block rotations, and thrust, strike-slip, and normal 
faulting have all been related to the convergence between two 
lithospheric plates having large continental masses. Much of the 
research on continental collisions has focused on the Indi a-Eurasia 
collision; however there are several other examples on continental 
collision which are ongoing today or have occurred in the past. 
Several models have been proposed to account for some or all of the 
features of continental collision. Frequently one of the colliding 
continents is more deformable than the other. A viscous sheet/ 
rigid punch model has enjoyed considerable success in accounting 
for many of the tectonic features for this case. For example, 
Tapponnier and Molnar (1977) applied this model in Asia to explain 
the extensive left lateral strike-slip motion that results in an 
eastward extrusion of a portion of China away from the collisional 
front. The model is also consistent with the right lateral fault 
motion observed in Southeast Asia, the compression observed north 
of Himalayas, the crustal thickening in the Tibet, and perhaps even 
the extensional deformation observed in Shansi graben and Baikal 
rift zones. England and McKenzie (1982) and Vilotte et al . (1982) 
have derived mathematical representations of this model and have 
done extensive numerical experiments to elicit the dependence of 
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the deformation on model parameters. Their results and those of 
Cohen and Morgan (1986) show that the critical model parameters are 
the stress exponent and multiplicative constant of the stress- 
strain rate constitutive law, the parameters expressing the rela- 
tive magnitude of the horizontal tectonic and i sostatical ly induced 
stresses, and the boundary conditions placed on the plate edges. In 
this study a finite element realization of the England and McKenzie 
(1982) model was used to study the modes of compression predicted 
to occur as a consequence of continental collision. 


RECENT ACCOMPLISHMENTS AND SIGNIFICANCE 

Figure 1 shows the essential features of the viscous sheet-rigid 
punch model . The lithosphere of the deformable continental plate 
is represented by a thin, incompressible viscous sheet governed by 
a nonlinear power law constitutive relation: 

r^ = beW"- 1 ^ (1) 

where r and e are the deviatoric stress and strain rate tensors, E 
is the second invariant of the strain rate tensor, and B and n are 
numerical parameters. A single constitutive law is used for the 
entire lithosphere but the division of the lithosphere into crust 
and mantle layers is accounted for in vertically averaging the 
deformation. During the collision a portion of the southern 
boundary of the deformable sheet is moved northward at a prescribed 
rate in order to represent the advance of a rigid punch or nonde- 
forming continental plate. As a consequence both horizontal defor- 
mations and crustal and mantle thickness changes occur in the 
sheet. Local isostasy is maintained by assuming Airy compensation 
against a crustless oceanic plate. Conservation of crustal mass is 
also maintained through an incompressibility condition. Figure 2 
shows an example of the deformation of a finite element mesh that 
is used to model the lithosphere. Symmetry is assumed about the 
Y-axis so that only the right-hand side of the deforming sheet 
needs to be shown. 

The major feature examined in this report is the contraction 
rate observed in front of the punch; particular emphasis is placed 
on its spatial and temporal dependence. Specifically e yy is 
plotted as a function of distance from the punch at various times 
along the axis X=0. A moving coordinate system is used where the 
value Y=0 is located on the boundary between the advancing punch 
and deforming plate. There are two classes of patterns for e y y. 
When the horizontal stresses are sufficiently large that they domi- 
nate the isostatic relaxation, the deformation pattern shows very 
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little time-dependence. For a linear rheology, the peak value of 
the strain rate occurs not on the punch boundary but at a consider- 
able distance into the plate interior. Conversely when the B value 
of the constitutive law is sufficiently small that the tectonic and 
isostatic stresses can be of comparable magnitude, the strain rate 
field can propagate ahead of the advancing punch. This propagation 
is more efficient for a nonlinear rheology with typical values 
being few cm/yr. Some specific examples of these two types of 
behavior are presented below. 

Figure 3 shows the dependence of the strain rate on distance from 
the punch for two difference times for a linear rheology with B 
equal to 1. x lO 2 ^ Pa sec (corresponding to a viscosity of 5. x 
10 23 Pa sec). The peak value of the strain rate occurs not at the 
punch boundary, but at about 1000 km from the punch, a value 
determined primarily by the punch width. The strain rate field is 
very broad. The value of the strain rate at the northern boundary 
of the sheet (Y=3600 km) at t=0 is still about 50 percent of the 
peak value. The broadness of the strain rate field for the linear 
rheology makes the results somewhat sensitive to the choice of 
distant boundary conditions (for these results, the northern boun- 
dary of the sheet is held fixed). Additionally, the strain rate 
field shows only minor changes with time. This general behavior 
persists with decreasing viscosity until the viscosity is smaller 
than that suggested by various geophysical measurements. However 
with a viscosity of 5. x 10 23 Pa sec the differences between the 
strain rates at the two times shown on Figure 4 are more 
substantial. Following the buildup of topography and crustal 
thickening, isostatic effects reduce the strain rate near the 
colli si onal boundary and increase them farther away. Still, the 
position of the peak strain rate advances only very slowly away 
from its initial position relative to the boundary. The results 
for a nonlinear rheologies are quite different. In Figure 5, for 
example, the values n=3 and B=5 x 10l 3 (cgs units) were used which 
correspond to the experimentally-determined values for typical 
lithospheric rocks. With the nonlinear rheology, the effective 
viscosity can be very small in the contact region and larger at a 
distance. As a consequence when the initial contact occurs between 
the plates, the strain rate is a maximum closer to the boundary 
than for the linear rheology. As the crust thickens with time, the 
stresses associated with the isostatic process increase and become 
comparable to the stresses due to the collision. When this condi- 
tion is achieved, further thickening of the region is resisted. 
Instead the crustal material is extruded northward (and eastward) 
resulting in an increase in the strain rate and crustal thickening 
in front of the punch. The position of the strain rate maximum 
propagates northward. The rate of propagation is the time deriva- 
tive of the extremum position. For the parameters used in this 
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calculation, the propagation velocity reaches approximately 3 
cm/yr. The calculation can be repeated for a nonlinear rheology 
with a higher B value or with no density contrast between the 
mantle and crust. These results are shown in Figure 6 which 
illustrates a case where the propagation effect has been eliminated 
with this change in parameter values. Two general conclusions can 
be reached by comparing these results. First, the propagation of 
the strain field is much more pronounced for n=3 than n=l, and 
second, the propagation results from an interaction of the tectonic 
and isostatic stresses and disappears when the gravitational 
effects vanish. 


FUTURE EMPHASIS 

The deformation features associated with continental collision are 
influenced both by the compressive and shear stresses generated by 
the collision and by the extensional stresses generated by the 
induced isostatic response. A fully satisfactory treatment of this 
problem requires a three-dimensional analysis where stress and 
strain rate can vary with depth. Future work will be directed 
toward the development of a complete three-dimensional finite 
element model capable of analyzing both the horizontal and 
vertical deformations. 
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Viscous sheet/rigid punch model of continental deformation. The sheet is an incompressible viscous fluid governed by a single 
linear or nonlinear constitutive relationship but with crustal and mantle material. A thin layer approximation is used with depth 
averaging of deformation and stress. In the finite element realization of this model the effects of the rigid punch are modeled by 
imposing northward motion on a portion of the southern boundary. 
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Figure 2. Deformation of the finite element mesh 24 million years subsequent to collision 
initiation. The grid contains nodes initially separated by 200 km. The punch half 
width is 1200 km; its northward velocity is 5 cm/yr. The numbers are crustal 
thicknesses in kilometers, with the initial crustal thickness equal to 35 km. The 
rheological law is the same as used in Figure 3a. 
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Figure 3. Strain rate versus distance from punch at collision initiation and after 24 million 
years. Linear rheology (n = 1). 

a) B = 1. x 10 24 Pa sec. 

b) B = 1. x 10 22 Pa sec. 
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Figure 3. Strain rate versus distance from punch at collision initiation and after 24 
million years. Linear rheology (n=l). B=l. x 10 24 Pa sec. 
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Figure 4. Strain rate versus distance from punch at collision initiation and after 24 
million years. Linear rheology (n=l). B=l. x 10 2 Pa sec. 
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Figure 6. Strain rate versus distance from punch at collision initiation and after 24 
million years. Nonlinear rheology (n=3). B is very large or crustal and 
mantle densities are equal. 



ANALYSIS OF INTRAPLATE DEFORMATION IN THE CENTRAL INDIAN BASIN 


Maria T. Zuber 


OBJECTIVE 

Congressional intraplate deformation of oceanic lithosphere is 
more clearly developed in the Central Indian Basin than anywhere 
else on the seafloor. This deformation violates the fundamental 
assumption in plate tectonics that lithospheric plates behave 
rigidly, deforming only along their boundaries. As most global 
dynamic and kinematic tectonic models are formulated on the basis 
of the rigid plate hypothesis, it is of obvious importance to 
characterize and understand those areas which are deforming 
internally. The purpose of this work is to investigate the nature 
of intraplate deformation of oceanic lithosphere using models of a 
compressing density and strength (viscosity) stratified medium. 
Model results are compared to the observed geophysical and structu- 
ral signature of the intraplate deformation in order to constrain 
the level of stress and style of sub-surface deformation in this 
part of the Indo-Australian plate. 


BACKGROUND 

Intraplate deformation in the Central Indian Basin, in the vicinity 
of the Ninetyeast Ridge, consists of approximately east-west tren- 
ding topographic undulations with an amplitude of 1-3 km and with a 
regular spacing of about 200 km (Weissel et al., 1980). The basement 
undulations are oriented approximately normal to the direction of 
maximum compressive stress in the Indian plate, which suggests that 
the deformation is related to forces resulting from the collision 
of India and Eurasia to form the Himalayas. Associated with the 
topography are large gravity and geoid anomalies, extensive 
seismicity, and high heat flow. A map illustrating the east-west 
trending anomalies in the marine geoid as derived from the Seasat 
and GEOS-3 altimeter data is shown in Figure 1. 

To better understand the nature of the intraplate deformation, 
models such as shown in Figure 2 have been developed (Zuber et al., 
1985a, b) in which the oceanic lithosphere is treated as a strong 
viscous or plastic layer of uniform strength and thickness h 
overlying a viscous substrate in which strength decreases with 
depth. This simple rheological representation, which approximates a 


37 


lithosphere in which deformation is accommodated by brittle frac- 
ture near the surface and thermally-activated ductile flow at 
greater depths, has been successfully applied to describe length 
scales of continental tectonic deformation (e.g., Fletcher and 
Hal let, 1983). The model lithosphere is assumed to be in a state 
of uniform horizontal compression. 

Two mechanisms are considered to explain the regular spacing of the 
seafloor topography, namely, flexural buckling of a thin plate and 
the hydrodynamic growth of compressional instabilities in a layered 
medium. In both models periodic deformation develops at a dominant 
wavelength ( ) , defined by the wave number k ( j(=2ir/X c i ) at which 
viscous resistance to deformation is minimized. The dominant 
wavelength is a function of the physical properties of the model 
lithosphere, in particular the layer thickness and the density and 
strength stratification of the medium. The dominant wavelength is 
determined from the solution for the instantaneous flow in the 
compressing layer and substrate. The model results, combined with 
the observed wavelength of deformation in the Indian Ocean, other 
geophysical observations, and experimental results on rock rheology, 
provide constraints on the mechanical and compositional structure 
of the oceanic lithosphere. 


RECENT ACCOMPLISHMENTS 

For a range of possible models, three modes of intraplate deforma- 
tion are possible. These are illustrated in Figure 3 for represen- 
tative values of the dimensionless parameter S, which relates the 
buoyancy force due to topographic variations to the vertically 
averaged strength of the layer. For S=0.1 deformation occurs by 
flexural folding, in which the strong layer of the lithosphere and 
the crust (which composes a fraction of the strong layer) retain 
approximately constant relative thickness. For S=1 deformation 
occurs by folding accompanied by thickening of the crust and layer 
in the vicinity of topographic highs. For S=10 deformation occurs 
in a pinch and swell mode. In the hydrodynamic instability model 
for a viscous layer, all three modes of deformation are possible, 
while if the layer is plastic, only the pinch and swell mode can 
occur. For the viscous buckling model, only the flexural folding 
mode of deformation is allowed. A comparison of the layer 
strengths and strong layer thicknesses associated with each model 
to the depths of intraplate seismicity and experimental results for 
rock rheology indicate that the S=0.1 and S=1.0 modes of deforma- 
tion are possible in the Central Indian Basin. Thus, on the basis 
of present data, either the hydrodynamic growth of instabilities or 
flexural buckling of a viscous layer may explain the internal 
deformation of the Indian plate. 
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FUTURE EMPHASIS 


To distinguish between the flexural folding (S=0.1) and layer 
thickening (S=l) modes of deformation requires measurements of 
crustal thickness associated with topographic highs and lows within 
the region of intraplate deformation. Seismic refraction profiling, 
which will allow the determination of crustal thickness, as well as 
the collection of more gravity and heat flow data in the Central 
Indian Basin are planned for the coming year. 
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OF POCnR QUALITY 



Figure 1. Shaded relief map showing east-west trending geoid anomalies correlated with seafloor 
deformation. The mean sea surface was derived from Seasat and GEOS-3 altimeter 
data (Marsh et al., 1984). 
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Figure 2. Model of compressing oceanic lithosphere with strong viscous or plastic layer of 

thickness h overlying weaker viscous substrate with decay depth f . The strong layer 
is composed of the crust and that part of the upper mantle which makes the greatest 
contribution to lithosphere strength. e xx is the mean horizontal strain rate, and p, r, 
and n are the density, strength, and stress exponent in the steady state flow law 
relating stress and strain rate. 
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strong layer 


Figure 3. Instantaneous displacement fields illustrating three possible modes of intraplate 
deformation. For deformation by the hydrodynamic growth of compressional 
instabilities all three modes are possible, while for deformation by viscous buckling of 
a thin plate only the S = 0.1 mode is allowed. The S = 0.1 and =1 modes are 
consistent with Indian Ocean intraplate seismicity and experimental rheological results. 
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CRUSTAL DEFORMATION UNDER THE TIBETAN PLATEAU AND THE 


HIMALAYAN RANGE 


Han-Shou Liu 


OBJECTIVE 

The purpose of this paper is to calculate the crustal stresses 
under the Tibetan Plateau and the Himalayan Range from plate tec- 
tonics and satellite-derived gravity data. Our intent is to use 
this crustal stress model as a geodynamical basis for interpreting 
seismotectonic characteristics and integrating the satellite- 
derived gravity data into the framework of plate tectonics. 


BACKGROUND 

Since Argand (1924) presented his ideas on the tectonics of Asia, 
Tibet has emerged as a symbol of continental collision. Today, 
within the framework of plate tectonics, there is little doubt that 
the massive plateau which maintains an elevation of 4.5 km, over 
600 x 1000 km^ north of the Himalayas, is anything but a direct 
consequence of the most spectacular India-Eurasia plate collision 
which has resulted from global plate motions. The forces that have 
pushed up this immense region must produce a crustal stress field 
in Tibet. Because a number of tectonic plates on the Earth move in 
a rather complicated way, the critical problem is how to derive the 
crustal stress field under Tibet from the stress functions of the 
global plate boundaries. In this paper mathematical methods for 
calculation of the stress field in the crust caused by global plate 
motions are developed. 

The methods of stress calculation developed in this paper are based 
on three fundamental premises: (1) Kelvin's theory for the 
deformation of spherical shells is correct; (2) Runcorn's theory 
for the geoid and mantle convection is correct; and (3) The rates of 
crustal creation and destruction along the global plate boundaries 
and the satellite-derived gravity data are correct. Given these 
premises, this study relates the geophysical theories to satellite 
observations. 


Despite decades of research by geoscientists, the crustal deforma- 
tion of the Tibetan Plateau and the Himalayan Range remains 
enigmatic. In this paper, a relatively underexploited approach to 
the study of stress fields which lead to crustal deformation is 
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presented. The exploration of the stress environment in the crust 
inferred from plate tectonics and satellite-derived gravity data 
would provide geophysical explanations for geological field 
observations and seismological dislocations in Tibet (Liu, 1985). 


RECENT ACCOMPLISHMENTS 

This paper uses plate tectonics and satellite-derived gravity data 
to further discussion of crustal deformation under the Tibetan 
Plateau. The first of three contributions is a spherical harmonic 
analysis of the global plate boundary system. A distribution of 
470 Dirac delta functions is applied to describe the generating 
forces according to the rates of crustal creation and destruction on 
the plate boundaries. Analysis of the extensional and compres- 
sional forces in the spreading and subducting zones shows that the 
present global plate motion causes compressional stresses in the 
N-S direction under the Tibetan Plateau (Fig. 1 and 2). The second 
contribution is the calculation of the crustal stresses in Tibet as 
inferred from satellite gravity data. By applying solutions to the 
problem of the spherical shells, the satellite-determined stresses 
indicate that the upwelling mantle material under Tibet induces N-S 
and E-W extension (Fig. 3). Finally, a superimposed stress system 
is constructed. This stress system shows that the present crustal 
deformation in Tibet does not produce N-S shortening but generates 
E-W extension (Fig. 4). 

The results of this paper have provided geophysical explanations 
for geological field observations in Tibet and fault plane solutions 
of earthquakes in the Tibetan side of the India-Eurasia collision. 
The stress patterns reveal that the cold downwelling mantle convec- 
tion flow beneath southern Tibet pulls the Indian plate down but 
applies a bending moment on the end of the plate to uplift and 
support the mass of the Himalayas. 


SIGNIFICANCE 

(1) Interpretation of crustal deformation 

Crustal deformation is one of the most puzzling aspects of 
plate tectonics. The geophysics governing crustal deformation are 
not yet known. At present there are two theoretical points of 
view. One is that deformations may be described by the motion of a 
number of tectonic plates moving in a rather complicated way. The 
other is that intraplate deformations are caused by the convection- 
generated subcrustal stress systems. These descriptions are now 
sufficiently precise to allow analyses to be carried out with 


44 



mathematical rigor. Therefore, the relationship of crustal 
deformations to the combined stress field caused by the global 

plate motions and mantle convection should be clear. 

If the momentum involved in crustal deformation under the Tibetan 
Plateau can be neglected, all crustal deformation in this region 
must be caused by forces now acting. There are two possible forces 
which could in principle maintain the observed deformation: force 

due to plate motion and force on the base of the crust which is 
generated by mantle convection. The analysis of stresses caused by 
the motion of a number of tectonic plates moving in a rather 
complicated way shows that the Tibetan crust is under N-S 

compression. The calculation of the crustal stresses caused by the 
subcrustal stress field as inferred from the satellite-derived 

gravity data shows that this region is under N-S and E-W extension. 

Therefore, crustal deformation under the Tibetan Plateau could be 
formed by the superposition of these two force systems, the N-S 
compression associated with the global plate motions plus the N-S 
and E-W extension caused by the upwelling mantle convection cell. 
In this case, the E-W component of the tensional stresses from 
mantle convection would produce E-W crustal extension and the N-S 
component of the tensional stresses from mantle convection just 
cancels the N-S compression associated with plate motion. The 
combined force system shows that the present crustal deformation in 
Tibet does not produce N-S shortening but generates E-W extension 
(See Fig. 4). 

The validity of these theoretical results is subject to observa- 
tional tests. The most important observational tests which verify 
this deformation model are: (1) fault plane solutions from earth- 

quakes in Tibet and (2) geological results from field expeditions 
in Tibet. 

Fault plane solutions for all crustal and intermediate depths events 
beneath the highest parts of the Tibetan Plateau show combinations 
of normal and strike-slip faulting with tensional stresses oriented 
approximately in the E-W direction (Molnar and Tapponnier, 1978; 
Molnar and Chen, 1983). None of these fault plane solutions from 
earthquakes exhibit components of thrust faulting. Therefore, 
seismic data confirm that active tectonics in Tibet are dominated 
by E-W extension. 

Tibet has been much studied by geologists in China. Important 
results of their studies have been extended by the 1980 French- 
Chi nese joint study of the geological structure, formation and 
evolution fo the Earth crust and upper mantle of this region 
(Tapponnier et al., 1981a). The most surprising field observations 
were the comparatively small or no amount of N-S shortening and 
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deformation. A pure E-W extension in Tibet seems to reflect a 
particular state of stress in the crust and upper mantle, where the 
upwelling mantle material as inferred from the satellite-derived 
gravity data can be the main driving force. Microtectonic observa- 
tions by Tapponnier et al . (1981b) support this inference. 

(2) Uplift of the Himalayas 

Figure 3 shows an upwelling mantle convection cell of lenticu- 
lar shape under the Tibetan Plateau. The general geological set- 
ting and continental subduction along the Himalayan Range are 
outlined in Figs. 5 and 6. Figures 5 and 6 show that the Main 
Central Thrust (MCT), the Main Boundary Thrust (MBT) and the 
Boundary Thrust Fault (BTF) (Gansser, 1964; Seeber et al., 1981) 
are in the compressional stress zone generated by mantle convection. 

By assuming that the topography of the Himalaya is supported 
by the Indian plate, Lyon-Caen and Molnar (1983) discussed the 
subduction of the Indian plate beneath the Himalayan Range. They 
have concluded that regardless of flexural rigidity of the plate, 
the weight of the material in the Greater Himalaya depresses the 
Moho too much to fit the observed gravity anomalies unless a large 
bending moment is applied to the north end of the plate. The origin 
of such a bending moment is unknown. 

The possible source of a large bending moment may be the torque 
applied to the Indian plate by the downwelling mantle flow which is 
relatively cold and dense under the range. If the under-thrusted 
Indian plate extends north of the Indus suture and is pushed for 
bending by the convective downwelling of the lower lithosphere, the 
force acting on the north end of the plate will apply a torque to the 
portion of the plate beneath the Himalaya. 

A cross section PT of the mantle convection pattern beneath 
the Himalaya as inferred from satellite gravity data is shown in 
Fig. 7. Figure 7 shows that there is cold and dense downwelling 
material beneath southern Tibet that pulls the Indian plate down 
slightly, but exerts a torque to the end of the plate so as to 
apply a bending moment. Therefore, the mantle convection cell of a 
lenticular shape under the Tibetan Plateau applies a bending moment 
on the end of the Indian plate to uplift and support the mass of 
the Himalayas. 


46 


FUTURE EMPHASIS 


The superimposed stress system based on plate tectonics and the 
satellite-derived gravity data seems to provide a geophysical basis 
for iterpreting the complicated deformation of the Himalayan Range 
and the Tibetan Plateau. Future work will concentrate on identi- 
fying more planetary tectonic stress regimes such as the Ishtar 
Terra on Venus and the Tharsis Plateau on Mars. 
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Figure 1. Distribution of .470 forcing functions along plate boundaries. + and • represent 

points of positive and negative delta functions in the spreading and subducting zones 
respectively. 
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generated by mantle convection. 
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Figure 7. Deformation of the Himalayas: Constraints from space geodynamics. 


CHAPTER II 


SEA SURFACE TOPOGRAPHY AND OCEAN/SOLID-EARTH INTERACTION 


OVERVIEW 

The oceans have a fundamental role in regulating major processes 
occurring on the earth's surface. Weather and climate are 
profoundly influenced by the heat energy stored and transported by 
the oceans, and the oceans are the primary source of water that 
reaches the continents as rain and snow. The immense importance of 
ocean transport is self-evident, and timely information and predic- 
tions of ocean surface conditions and currents are also important 
to both marine shipping and population concentrations along coasts. 
A basic element in understanding these practical aspects of oceano- 
graphy is ocean dynamics which embraces the general ocean circula- 
tion, tides, current systems, the air-sea interface, sea state, 
storm surges, etc. An essential source of information is the 
topography of the ocean surface. Ocean currents have a topographic 
signature — a slope of the sea surface perpendicular to the flow 
caused by the interaction of the current with the Coriolis force. 
For example, the Gulf Stream has a rise of approximately 1 meter 
over a distance 50 km. Satellite altimeters measure the distance 
between the satellite and the ocean surface and thus can obtain 
information on the topography of the ocean surface along the sub- 
satellite track. Recovery of the sea-surface topography requires 
knowledge of the satellite orbit. In order to recover oceanographic 
and meteorological information, it is necessary to separate the 
mean sea surface topography (geoid undulations) from the measured 
topography obtained with satellite altimetry. The overall objec- 
tive of the research described in this chapter is to apply remote 
sensing data exemplified by satellite altimetry, satellite laser 
ranging, orbit perturbation, and gravity data as the basic data 
types used to perform studies of dynamic ocean and earth processes. 
Research programs are underway in the following areas, to analyze 
and interpret altimetry data, spaceborne radiometer data, and 
satellite laser ranging data from Lageos: (1) The computation of a 
global mean sea surface based upon a combination of GE0S-3 and 
Seasat altimeter data. (2) Examination of large scale circulation 
of the world's oceans utilizing remotely-sensed satellite-borne 
radiometer in conjunction with in situ data. (3) Development of 
mathematical models and computational algorithms for the computa- 
tion of any tide constituent and its esimated error from the major 
tides and their associated errors. (4) Develop and test an inter- 
pretation technique which allows accurate extrapolation of tidal 
height fields in the ocean basins. (5) Develop numerical modeling 
to examine the effect of the wind-driven oceans on the rotation of 
the earth. (6) Study the effects of the 13-month pole tide on the 
earth's rotation. (7) Review and validate the Seasat altimeter 
height calibration technique. 

Contributors to this chapter are: Demosthenes C. Christodoul idi s , 
Ronald Kolenkiewicz, Chester J. Koblinsky, James G. Marsh, William 
P. O'Connor, Braulio V. Sanchez and Daniel J. Steinberg. 
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GLOBAL MEAN SEA SURFACE BASED UPON 
A COMBINATION OF GEOS-3 AND SEASAT ALTIMETER DATA 


James G. Marsh 


OBJECTIVE 

The objective of this work has been the computation of a global 
mean sea surface map based upon the GEOS-3 and SEASAT altimeter 
data sets. The dominant error source in the satellite altimeter 
data is the radial orbit error. This error is primarily due to 
earth gravity model errors. These errors are long wavelength in 
nature, primarily once per revolution. In order to minimize the 
effects of the radial orbit error, a combination of crossing-arc 
techniques and accurate reference orbits have been used. 


BACKGROUND 

The mean sea surface provides a reference surface for the detection 
of mesoscale ocean circulation features, global ocean circulation 
patterns and detailed information over the oceans on the internal 
structure of the earth. 

The GEOS-3 and SEASAT altimeter experiments have provided important 
sets of data for the computation of a mean sea surface. The 
precision of the GEOS-3 data is about 30 cm and the precision of 
the SEASAT data is 5-10 cm for length scales of a few thousand 
kilometers. The dominant frequency of the gravity model errors in 
the radial position of the satellite is once per orbital revolution 
and can thus be wel 1 -represented by a linear function of time over 
distances of a few thousand kilometers. In the present analyses 
crossing-arc adjustment techniques were used to solve for linear 
trend parameters in order to absorb radial orbit errors. The linear 
trend parameters also tend to absorb other errors, e.g., tropos- 
pheric and ionospheric refraction errors, and tidal model errors. 
Any geographically correlated gravity model errors cannot be 
removed by this procedure and thus absolute errors of a few deci- 
meters over ocean basin scales may still exist in the mean sea 
surface. 


PRECEDING 
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RECENT ACCOMPLISHMENTS 


Previous mean sea surface maps had been computed at GSFC on grid 
spacings of 0.25°x0.25° and 0.5°x0.5°. Detailed analyses indicated 
that the along-track data density and the cross-track spacing 
supported a further reduction in the grid spacing. In the present 
analyses a grid spacing of 1/8° was used. The limiting factors in 
the computations are now due to the mesoscale ocean height varia- 
bility due to the major current systems. In the oceanic quiet 
areas, the rms crossovers were usually less than 10 cm, while over 
the major Western Boundary Current systems residuals of 20 cm to 30 
cm were observed. 

The data analysis technique consisted of: the computation of 65 

regional crossover solutions which covered the global set of 
altimeter data; the adjustment of these regional solutions into a 
global reference grid, and the subsequent adjustment of the GEOS-3 
data set on a pass-by-pass basis into this SEASAT reference grid. 
A bi-quadratic local surface fitting function was used to determine 
a sea height value representative of the sea surface heights at 
each grid location. A grid spacing of 0.125° was used for the 
computations . 

This gridded data set has been analyzed using image processing 
techniques to enhance the short wavelength topographic features 
which are present in the map. The map is presented in the figure 
at the end of this section. Such features as the mid-Atlantic 
Ridge and the detailed fracture zone system, fracture zones in the 
Pacific (e.g., Mendocino, Murray and Clipperton and others), sea 
mount chains and other bathymetric features are clearly shown. In 
order to emphasize the short wavelength features, a long wavelength 
geoid corresponding to the GEM 10B (12,12) model has been removed 
from the mean sea surface. 


SIGNIFICANCE 

The shaded relief maps such as the one shown in the figure provide 
an important new basis for the definition of ocean circulation 
processes. In addition, an important basis for the interpretation 
of geophysical features over the oceans has been provided. 


FUTURE EMPHASIS 

Significant progress is being made in the refinement of the geodetic 
models for the earth gravity model, the reference coordinate system, 
earth and ocean tides and other parameters. These new models are 
expected to provide a major reduction in the radial orbit errors 
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for SEASAT. These new models will be used in combination with the 
SEASAT Doppler tracking data in order to form a more accurate 
reference network for the adjustment of the altimeter data. 

Negotiations are underway for the acquisition of the GEOSAT alti- 
meter and tracking data. The higher quality of these data and 

more favorable area to mass ratio versus SEASAT are expected to 
provide a significant improvement in the overall accuracy of the 
mean sea surface. It is anticipated that the amplitude of the 

geographically correlated gravity model errors can be reduced 
when the new TOPEX gravity model is used. 
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Figure 1 . Mean sea surface based upon SEASAT and GEOS-3 satellite altimeter data 


OCEAN CIRCULATION STUDIES 


Chester J. Kobl insky 


OBJECTIVE 

The purpose of these studies is to utilize remotely-sensed signa- 
tures of ocean surface characteristics from active and passive 
satelliteborne radiometers in conjunction with in situ data to 
examine the large scale low frequency circulation of the world's 
oceans . 


BACKGROUND 

For the past several years oceanographers have been able to study 
the ocean using a variety of passive and active radars in space on 
such missions as GEOS, NIMBUS, TIROS, and SEASAT. It is now 
possible to usefully measure a number of sea surface variables from 
space including temperature, winds, and dynamic topography. Such 
measurements complement observations of the ocean interior to 
significantly improve the understanding of marine physics including 
the dynamics of mesoscale processes, the momentum exchange between 
atmosphere and ocean, and the general circulation. 

Mesoscale processes are ubiquitous in the ocean. They have length 
scales of approximately 10 to 1000 km and time scales of days to 
months. In the deep ocean, isolated vortices or eddies such as 
western boundary current rings are an important mesoscale process. 
These features can have a dynamic topography anomaly of 25 to 100 
dyn cm that is observable with a satellite altimeter and can align 
the sea surface temperature field such that the eddy can be 'seen' 
with a satellite radiometer. Eddies can rotate rapidly, up to one 
revolution in a few days, and move about, typically 5 kilometers a 
day so that large scale rapid sampling from satellites provides a 
unique synoptic view of these processes. 

Direct forcing of the ocean by the atmosphere through the stress of 
the wind on the sea surface provides an important exchange of global 
angular momentum. Direct observations of low frequency deep-ocean 
currents have now been made that demonstrate the ocean's response 
to wind forcing (Niiler and Koblinsky, 1985). Low frequency wind- 
driven deep-ocean currents are thought to provide a background 
level of eddy kinetic energy and drive the general circulation. 
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On climatological scales the general circulation of the ocean and 
its variability is important to an understanding of the global heat 
and momentum balance and biogeochemical processes. The capability 
of satellite altimetry to observe the global sea surface dynamic 
topography now allows one to significantly improve the kowledge of 
the general circulation. Furthermore, the large scale ocean 
circulation must provide a couple between the atmosphere and solid- 
earth that may affect polar motion, earth rotation, and satellite 
orbits. 


RECENT ACCOMPLISHMENTS AND SIGNIFICANCE 

Two studies have been submitted which examine mesoscale processes 
in the California Current System with in situ and satellite 
observations (Simpson, et al., 1986; Haury, et al., 1986). These 
papers utilize a long time series (3 years) of west coast imagery 
from the NIMBUS-7 Coastal Zone Color Scanner (CZCS) and the NOAA 
Advanced Very High Resolution Radiometer (AVHRR) to study time 
changes of sea surface pigment and temperature. Over 30 years of 
hydrographic and biological samples from the California Cooperative 
Fisheries Investigations were used to examine long term change of 
physical, chemical, and biological in the upper 500 meters of the 
ocean and relate them to surface changes from the satellite imagery. 
In addition, comparisons were made with a very intensive ship survey 
made in January 1981. It is demonstrated that the offshore eddies 
in the California Current play a major role in the offshore movement 
and mixing of coastal waters. 

Work is continuing on studies of the wind-forced circulation of the 
North Pacific. One study is exploring the geographic variability 
of the mesoscale wind-forced response of the ocean from over 30 
current meter moorings, each extending from 1 to 3 years in duration 
(Koblinsky, Niiler, and Schmitz, 1986). This study clearly shows a 
strong ocean response in the winter months because of the intensifi- 
cation of the North Pacific winter winds. A second study with J. 
Price (WHO I ) and P. Niiler (SIO) is examining the generation of 
inertial/internal waves by large scale, intense storms that frequent 
the eastern North Pacific in the winter. 


FUTURE EMPHASIS 
A. Mesoscale Processes 

Research will continue on the merger of several different ocean 
remote sensing measurements with shipboard observations to under- 
stand ocean mesoscale processes. 

1) An effort to understand the quantitative description of mesoscale 
ocean eddy fields by satellite altimeters has been initiated using 
eddy resolving ocean numerical models. 
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2) J. Marsh (Geodynamics Branch), F, Sci remammano (RIT) and the 
author are continuing to study mesoscale variability in the SEASAT 
and GEOS-3 altimeter data sets. Techniques have been developed by 
Sci remammano to extract mesoscale features according to wave number. 
Comparative studies of individual features in eddy intense regions 
of the ocean are in progress. Techniques for extracting larger 
scale variability are being pursued using the long time series of 
altimetry data in the western North Atlantic. 

B. Atmospheric-Ocean-Solid Earth Momentum Exchange 

1) A joint research project is being coordinated in an attempt to 
measure the mesoscale wind-forced response of the North Pacific 
using satellite altimeter data from the GEOSAT mission. These 
measurements will be compared/assimilated with in situ measurements 
of bottom pressure, currents and accoustic tomography being carried 
out by investigators at Scripps Institution of Oceanography. 

2) The numerical ocean general circulation model of Cox and Bryan 
has been obtained from NOAA/GFOL and implemented on the GSFC Cyber 
205 by Dan Steinberg, NASA Graduate Fellow. This model is being 
used to study the wind-forced circulation of the global ocean using 
estimates of the wind field over the past 10 years obtained from 
the National Meteorological Center. 

3) Maps of the global intensity and variation of mesoscale wind- 
forced ocean currents using a model function derived from direct 
observations and winds from the SEASAT scatterometer and the 
National Meteorological Center are being produced. 
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PREDICTION OF MINOR OCEAN TIDE CONSTITUENTS 


Demosthenes C. Christodoul idis 


OBJECTIVE 

The most comprehensive models for ocean tides currently available 
represent only the major tides. Whereas these constituents account 
for over 90% of the global tidal amplitude, the remaining minor 
constituents can have significant perturbing effects on a 
satellite's orbit. The major tides have been determined by others 
through accurate numerical solutions of the Laplace Tide Equations, 
requiring a heavy computational burden. A simple predictive 
algorithm was developed for the computation of any tide constituent 
and its estimated error from the major tides and their associated 
errors. This process provides oceanographically reasonable 
a priori estimates of the minor tides and corresponding estimates 
of their errors which are necessary for the TOPEX gravity model 
improvement and for related precision orbit determination error 
analyses. 


BACKGROUND 

Modeling ocean tides has been a significant problem for precision 
satellite geodesy and geodynamics since the late 1960's. The 
effects on satellite orbits are often at the meter level and can 
reach tens-of-meters in satellite position. For the T0PEX/P0SEID0N 
mission, which is expected to be launched in February of 1991, a 
major goal is to map the global ocean topography to an accuracy of 
better than a decimeter — the total error budget is a decimeter in 
satellite radial position. The oceanic tides must be accounted 
for, both in their effect on the satellite orbit and in their effect 
on the altimetric measurement of the sea surface. 

Table 1 shows the estimated radial perturbation amplitude due to 
the major ocean tide constituents on the proposed TOPEX orbit and 
on the GEOS-3 orbit. This analysis was based upon a Kaula type 
first order linear orbit perturbation theory. More than half of 
the constituents have effects which exceed 1 decimeter radially. 
These terms must be modeled. It is probable that the associated 
minor tides for some of these also must be modeled if the minor 
tide response is proportional to the tide raising potential of the 
major tide. 


64 


The pervasiveness of the ocean tide effect on satellites is readily 
shown. A crude estimate of the effect of the ocean tide is about 
10% of the body tide. A set of 53 satellite orbits was evaluated 
for their nominal ocean tide perturbations of 230 tidal frequencies. 
These 53 orbits represent a variety of orbital inclinations and 
altitudes, and all have reasonable tracking data histories. Figure 
1 shows the number of satellites having effects over .001 arcsec in 
the inclination as a function of tidal frequency. Satellites were 
also included if the principal third degree terms from an ocean tide 
decomposition produced a perturbation in the orbit eccentricity 
greater than 1 ppm. In this analysis, the amplitude of the ocean 
tide coefficients was assumed to be 1 cm. Note that the criterion 
of 1 ppm perturbation in the eccentricity is equivalent to the cri- 
terion of a .001 arcsec perturbation in the inclination. 

Of the 230 tidal frequencies, significant perturbations can be 
anticipated for about 150 of them. Numerical tide models are 
available only for the major tides, shown by the shaded bars, but 
not for any of the sideband or minor tide constituents. Thus, 
there are possibly about 140 minor tide constituents for which 
models are required in precision computations. Because of possible 
commensurabilities between the orbital periodicities of a satellite 
and the periods of the tides at the surface of the Earth, these minor 
tides may actually produce the larger perturbation of the satellite 
orbit. 

The ocean tide effects are thus an important concern for precision 
orbit calculations and also for geodynamic parameter recovery from 
analyses of satellite tracking data. Neglect of these effects may 
have significant impact on gravity model determination because the 
time varying potential of the tides has a similar form to that of 
the long wavelength geopotential. The best gravity field model 
available prior to the T0PEX modeling efforts is estimated to 
produce radial rms errors on T0PEX of 54 cm (GEM-L2 evaluation by 
Rosborough, 1985). Among the known modeling deficiencies in the 
development of GEM-L2 is that no comprehensive tide modeling was 
available. T0PEX requires a major improvement in the knowledge of 
the geopotential. Proper modeling the ocean tides is anticipated 
to make a major contribution to this goal. 

Solutions of the Laplace Tidal Equations require a large computa- 
tional burden. Oceanographic numerical models such as those of 
Schwiderski or Parke exist for the major tide constituents only. 
It was therefore necessary to develop a method which would provide 
reasonable estimates for the minor tides and corresponding esti- 
mates of their errors from the available major tide constituent 
model s . 
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RECENT ACCOMPLISHMENTS 


A method has been developed for predicting minor ocean tide 
constituents. This method has been used to derive models and error 
estimates for some 36 minor tides. The method is based upon 
interpolating the admittance of each desired minor tide from the 
admittances of the major tides. The use of the admittance was 
motivated by the study of Munk and Cartwright (1966). The 
Schwiderski tide models from the Naval Surface Weaspons Center 
(NSWC), which are the most complete set available, were adopted for 
this effort. Minor constituent models are generated on a 1 degree 
global grid matching that of Schwiderski and have also been con- 
verted to spherical harmonics for subsequent satellite studies. 

The semidiurnal, diurnal, and long period tidal bands are sepa- 
rately analyzed so that the range of frequencies for interpolation 
is limited. The procedure assumes that the tidal response is 
locally linear, i.e., within each band at each latitude, longitude 
point on the Earth's surface. This linearity assumption was 
adopted because global nonlinearities are expected to be small, and 
also for the practical reason that there are at best four points to 
interpolate over (or extrapolate from) in each tidal band. The fit 
is weighted according to the nominal errors in each tide at the 
given location. The nominal errors in the diurnal and semidiurnal 
band were global values supplied by Schwiderski. Nominal errors 
for the long period band were estimated as being proportionally as 
well determined relative to the equilibrium tide as M2, i.e. 12.8%. 
Note that the NSWC Mm and Mf tide amplitudes are smaller by a 
factor of 3 or 4 compared to their nominal equilibrium values. 
This suggests a conflict with the assumption of linearity of 
response across the long period tidal band. Proportionally, there 
is a much greater span of frequency variation in the long period 
band than in the diurnal or semidiurnal band. However, only three 
long period tides are availale so this frequency band cannot be 
further segmented to reduce the range of interpolation. 

Table 2 shows the global statistics associated with the local linear 
fits for each of the major tides. With the exception of 01 and M2, 
the rms global fit in each of the semidiurnal and diurnal bands 
shows that the linear model disagrees with the NSWC input by 
approximately the estimated error in NSWC. M2, the worst case, has 
a weighted rms disagreement of 7 cm out of a total 30 cm, which is 
an error in power of less than 5%. The fits in the long period 
band show that the long period tides are probably not adequately 
modeled with this procedure, in that the weighted rms residual 
amplitude is of the order of the entire NSWC rms tide amplitude. 
The standard deviations of unit weight given in Table 2 provide the 
factors by which the NSWC rms amplitude errors need to be adjusted 
in order to map the weighted residuals into the unit normal 
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distribution. The semidiurnal and diurnal bands are near unity, 
but the long period band is off by a factor of 2.5. Thus the linear 
model is not inconsistent with the semidiurnal and diurnal data, but 
it is inconsistent for the long period tides assuming the projected 
error estimates for these tides are correct. However, rms fits are 
still less than 40% of the equilibrium amplitude for these long 
period tides, and are only about twice the estimated nominal error 
in these tides. 

Figure 2 compares the global amplitude of the M2 tide of NSWC with 
that from the numerical method used. Figure 3 shows the amplitude 
of the total error (vector magnitude) as a function of latitude and 
longitude and also the corresponding percent relative error. The 
models are qualitatively the same, which indicates that this tide 
was not seriously corrupted by the procedure used. This is so for 
all of the semidiurnal and diurnal tides. Similar comparisons for 
the long period tides are less successful as is expected from the 
above discussion. These tides show regions of over 100% 
differences. 


FUTURE EMPHASIS 

An extensive comparison for each of the major tides in terms of 
amplitude, phase, admittance, and estimated errors is in 
preparation. Future emphasis will be in terms of using these 
models as error models for satellite precision orbit determination 
studies and for comparisons with spherical harmonic tide models 
recovered as part of the TOPEX gravity field improvement effort. 
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Table 1. Ocean Tide Radial Perturbation Amplitudes on Satellite Orbits 
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Figure 1. Histogram of the number of satellites with strong sensitivity to various tidal frequencies. The tidal frequency axis is not linear. The 
1 1 major constituent tidal families are marked by brackets. Numerical tide models exist only for the principal tidal frequency of 
each family shown by the black bars. 
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Figure 2. Comparison of M 2 Tide Amplitudes 
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Figure 3. Error in Interpolated M 2 Tide 
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AN OBJECTIVE ANALYSIS TECHNIQUE FOR EXTRAPOLATING TIDAL FIELDS 


Braulio V. Sanchez 


OBJECTIVE 

This work develops and tests an interpolation technique which 
allows accurate extrapolation of tidal height fields in the ocean 
basins by making use of selected satellite altimetry measurements 
and/or conventional gauge measurements. 


BACKGROUND 

The theoretical foundation of the method used in the calculations 
is based on Proudman's theory (1918) as reformulated by D.B. Rao 
(1966). The theory provides a formulation for calculation of 
gravitational normal modes and rotational normal modes of irregu- 
larly shaped basins with realistic bathymetry. The objective 
analysis uses orthogonal functions which form the basis to represent 
the water level fluctuation field. The' method has been applied 
successfully to compute the M2 tide in Lake Superior by Sanchez, 
Rao and Wolfson (1983). 


RECENT ACCOMPLISHMENTS 

A normal mode solution for the Atlantic and Indian Oceans has been 
obtained by means of 3°x3° finite difference grid. There are 1915 
velocity potential and 1717 stream function points in the grid. An 
eigenvalue solution was obtained for both fields. The normal mode 
solution was obtained by including 150 eigenfunctions from each 
field (velocity potential and stream function) into the Laplace 
tidal equations, a third eigenvalue solution is then obtained which 
yields the normal modes. The normal modes were used to obtain the 
coefficients for the forced solution and to compute power spectrums 
for the tidal components M2 and Kl. Table 1 below gives a list of 
the first 20 gravitational modes in order of decreasing period, for 
each mode the ratio of energies i$ given as well as the percentage 
of the total power which that mode contributes to the M2 and Kl 
tidal components. The three most powerful modes in the M2 spectrum 
are the ones with periods of 12.68 hours (7.73%), 11.10 hours (5.76%) 
and 11.55 hours (5.23%). Correspondingly the three most powerful 
modes in the Kl spectrum have periods of 22.89 hours (10.47%), 
19.13 hours (6.88%) and 20.61 hours (6.57%). 

The coefficients for the forced solutions yield the tidal amplitudes 
and phases when used in conjunction with the velocity potential 
eigenfunctions previously computed. These functions are linearly 
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combined to obtain the amplitude and phase. The computation of the 
forced solution includes the effects of the yielding of the solid 
earth to the tide generating forces through the Love numbers k 2 and 
h 2 ; the effects of self-attraction of the tide, tidal loading and 
bottom friction have not been incorporated in the theoretical 
simulation at this point. The amplitudes and phases have been 
computed at the velocity potential points defined in the 3°x3° 
grid. Figures 1 and 2 below give the structure for the M2 and K1 
tides in the Atlantic-Indian Oceans. The investigators involved 
are B.V. Sanchez, D.B. Rao and S. Steenrod. 


SIGNIFICANCE AND FUTURE EMPHASIS 

A method has been developed to objective analyze tidal amplitudes 
and phases to map the tidal constituents over an entire basin. The 
method allows the computation of the normal modes of the basin as 
well as a theoretical forced solution for each of the tidal 
components. The method has been applied to the Lake Superior basin 
(1983) and to the Atlantic-Indian and Pacific basins using a 6°x6° 
grid (1983,1984). 

The future emphasis will lie in the development of a 3°x3° solution 
in the Pacific Ocean and the application of the objective analysis 
using the 3°x3° solutions. 


REFERENCES 

Proudman, J., 1918. On the Dynamical Equations of the Tides, Proc. 
of the Lond. Math. Soc. 18, 1-68. 

Rao, D.B. , 1966. Free Gravitational Oscillations in Rotating 

Rectangular Basins, J. Fluid Mechanics 25, 523-555. 

Sanchez, B.V., D.B. Rao and P.G. Wolfson, 1983. Objective Analysis 
for Tides in a Closed Basin, Marine Geodesy , Vol . 9, NO. 1, 71-93, 
1985. 

Geodynamics Branch Annual Report - 1983, NASA TM 86123. 

Geodynamics Branch Annual Report - 1984, NASA TM 86223. 


74 



Table 1. Normal modes and power spectrum for the M2 and K1 components. Atlantic-Indian, 

3° x 3° grid. 


Period 

Energy Ratio 

Power % 

Power % 

(hours) 

K/P 

M2 

K1 

68.14 

1.66 

0.15 

2.53 

43.99 

0.93 

0.30 

1.27 

31.16 

0.79 

0.04 

1.81 

29.20 

0.94 

0.01 

2.06 

28.20 

0.94 

0.39 

2.09 

22.89 

0.90 

1.13 

10.47 

20.61 

0.87 

1.70 

6.57 

19.13 

1.20 

0.36 

6.88 

17.59 

0.99 

1.16 

4.43 

16.21 

0.92 

1.09 

4.41 

15.51 

1.05 

1.03 

1.42 

14.79 

0.89 

2.10 

0.76 

13.68 

0.98 

1.88 

4.50 

12.85 

1.05 

1.90 

4.12 

12.68 

1.05 

7.73 

1.41 

12.05 

1.15 

1.94 

2.42 

11.88 

1.03 

2.97 

0.88 

11.55 

1.03 

5.23 

2.01 

11.10 

1.01 

5.76 

1.12 

10.96 

1.11 

5.08 

2.47 
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OCEANIC EXCITATION OF THE EARTH'S ROTATION DURING EL NINO 


Daniel J. Steinberg 
Chester J. Kobl insky 


OBJECTIVE 

This work will utilize numerical modeling to examine the effect of 
the wind-driven oceans on the rotation of the earth. 


BACKGROUND 

The Chandler wobble is the earth's natural wobble. As a result of 
dynamic response of the oceans as well as the anelastic behavior of 
the mantle, the Chandler wobble should, eventually, be damped out. 
However, it has been observed continuously for many years, at times 
even increasing in magnitude rather than decreasing, indicating a 
recurrent source of excitation. The wobble can be excited by 
torques exerted on the solid earth externally or from deformations 
and motions within the earth. The most popular candidates for 
excitation are earthquakes and atmospheric forcing, but neither has 
been shown to have the necessary magnitude to match the observed 
excitation of the wobble. The oceans have been studied as a possible 
excitation source using simple numerical ocean models (Wahr, 1983). 
The authors propose to extend this research using a numerical 
general circulation model forced by observed global wind fields. 

Recent studies have shown a correlation between Southern 
Oscillation/El Nino phenomena and the polar motion (e.g., Chao, 
1985). The 1982-83 El Nino was the strongest observed in this 
century. It was responsible for large changes in the length-of-day 
(e.g., Rosen and Salstein, 1984). This correlation was due largely 
to the change in atmospheric angular momentum associated with the 
Southern Oscillation. The wind speeds and, therefore, the wind 
stress in the Pacific were enhanced during the 1982-83 El Nino 
(Wyrtki , 1984). This may have had a large effect on the ocean 
circulation. Such circulation must be modeled since data to estab- 
lish it are very sparse. This significant event in the short term 
climate will be the focus of this study. 


RECENT ACCOMPLISHMENTS AND SIGNIFICANCE 

The Cox and Bryan (1983) model has been obtained from NOAA/GFDL and 
implemented on the SESCC CYBER 205. Successful preliminary 
comparisons have been made between the Cox and Bryan model, a 
Sverdrup model and simple analytical models of wind drive circula- 
tion. Studies of variable wind forcing are now in progress using 
climatological wind fields. 
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FUTURE EMPHASIS 


Ocean models will be subject to the observed winds during the 1982- 
1983 El Nino period. The authors will be trying to determine how 
large the angular momentum and mass redistributions of a model 
ocean are when subject to realistic wind stress. In addition, how 
large an effect, such a wind-driven ocean can have on the length-of- 
day and on the Chandler wobble will be studied. The latter will be 
accomplished by a correlation between the excitation function 
derived from the model and an astronomical excitation function 
derived from LA6E0S polar motion data. 
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THE POLE TIDE AND THE EARTH'S ROTATION 


William P. O'Connor 


OBJECTIVE 

The objective of the research is to explain why the 14-month pole 
tide is enhanced above its equilibrium value in some coastal 
regions, especially along the Dutch and Danish coasts of the North 
Sea. 


BACKGROUND 

The 14-month Chandler wobble of the earth sets up the associated 
pole tide in the ocean, which has a maximum amplitude of a half 
centimeter in mid-latitudes. Theoretical work with Laplace's tidal 
equations has shown that the global structure of the pole tide 
should be very close to the equilibrium tide value determined by 
the wobble forcing potential, even when the effects of ocean 
boundaries are considered (O'Connor and Starr, 1983; O'Connor, 
1986). Since the amplitude of the pole tide is so small, with most 
of the data taken at coastal stations, it is difficult to determine 
the global structure of the pole tide from observations. The 
statistical analysis of long period sea level data records detects 
the pole tide signal, and shows that in some coastal regions the 
tide height is several times that value which would be predicted by 
equilibrium theory alone. This is especially true in the North and 
Baltic Seas. In particular there is a steep gradient of the observed 
pole tide height along the Dutch and Danish coasts, increasing 
northwards to a value five times that predicted by equilibrium 
theory. 

It is important to find an explanation for any significant deviation 
of the tide from equilibrium. One important question of geodynamics 
concerns the possibility of damping of the earth's Chandler wobble 
energy, and whether this could take place in the oceans. This 
damping can only occur when the pole tide departs significantly 
from equilibrium over a large part of the ocean basins (Lambeck, 
1980; Smith and Dahlen, 1981). Although these large deviations of 
the pole tide from equilibrium in coastal regions are not likely to 
result from large scale deviations of the tide from equilibrium 
over the ocean basins, it is certainly desirable to find an 
explanation for them. 
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The analysis done so far on the pole tide in the North Sea (Wunsch, 
1974) has shown that the slope of the bottom topography, deepening 
to the north, causes an eastward intensification of the flow regime 
and tide height. However, this effect alone is not sufficient to 
explain the large gradient in the tide height along the Danish 
coast. It is assumed that since the area of the North Sea is small 
compared to the global scale tidal potential, the direct effect of 
the tidal forcing on this sea is negligible, and so the Laplace's 
tidal equations contain no tidal forcing terms. Then the tidal 
forcing in the North Sea appears only as the boundary condition at 
the northern opening to the deep ocean. This has been shown to be 
an acceptable assumption in modeling the effects of the diurnal 
tides in the North Sea. 


ACCOMPLISHMENTS AND SIGNIFICANCE 

The work now in progress departs from previous works in that other 
non-tidal forcing mechanisms are considered. It is known that the 
North Atlantic and North European regions have a 14-month oscilla- 
tion in sea level pressure on the order of a few tenths of a 
millibar. These oscillations may or may not be an effect of the 
earth's Chandler wobble on the atmosphere-ocean system, but what- 
ever their origin, they have been detected in data records. These 
long period pressure oscillations must in turn give rise to low 
level atmospheric motion through the geostropic relationship giving 
a balance between the horizontal pressure gradient and coriolis 
forces. Then the wind stress at the air-ocean interface transfers 
horizontal momentum from the atmosphere to the ocean, and forces 
motion at least in the uppermost layers of the ocean. The effect 
of an explicit wind stress forcing term is included in the 
Laplace's equations of motion for a North Sea model. 

The wind stress forcing pattern at this low frequency must be 
obtained. An indirect approach is to use the atmospheric pressure 
data showing this small oscillation at the Chandler frequency 
(Bryson and Starr, 1977), and then compute the wind stress by the 
geostrophic relationship and momentum transfer equations. Histori- 
cally, most of the published wind stress data over the oceans has 
been given as monthly means over the course of an average annual 
cycle. Recently, however, Hellerman and Rosenstein (1983) have 
compiled a time series of monthly mean wind stresses over the 
oceans, and the analysis in progress seeks to determine directly 
the 14-month oscillation in wind stress over the North Sea. 
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The atmospheric momentum transferred to the oceans is proportional 
to the square of the wind speed. Then the oscillation in momentum 
transfer from air to ocean is proportional to the difference of the 
squares (as opposed to the square of the difference) of the 
oscillating wind speed from the average. Thus a small wind stress 
oscillation can transfer significant momentum to the sea, and drive 
a circultion in a shallow sea. It is this circulation that could 
cause the observed deviations of the pole tide from equilibrium. 


FUTURE EMPHASIS 

Work will continue on the low frequency dynamics of the ocean and 
the earth's rotation. This includes further work on low frequency 
meteorological forcing of the oceans. When these effects on sea 
level are better understood, it may be possible to determine whether 
the pole tide is in equilibrium over most of the oceans. 
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REEVALUATION OF THE SEASAT ALTIMETER HEIGHT CALIBRATION 


Ronald Kolenkiewicz 


OBJECTIVE 

The objective of this work is to review the SEASAT altimeter height 
calibration performed by NASA/GSFC (Kolenkiewicz and Martin, 1982) 
and to verify that the results are valid. These results are to be 
conveyed to NASA Headquarters Oceans Branch, the SEASAT and the 
TOPEX Project Offices at Jet Propulsion Lab (JPL), Naval Surface 
Weapons Center (NSWC) Dahlgren, VA, and Defense Mapping Agency 
Aerospace Center ( DMAAC ) . 


BACKGROUND 

During the summer of 1984 DMAAC (based on the analysis of their 
new World Geodetic System (WGS 84) results) concluded that the NSWC 
SEASAT-derived geoid was in error by about 2 meters. In an attempt 
to explain this discrepancy, NSWC found three potential error 
sources (Duke 1984): 

1. A NSWC software error — amounting to 71 to 80 cm. 

2. Neglect of the altimeter feedhorn-antenna separation distance 
in calculating the altimeter correction — amounting to 74 cm. 

3. A speculation of aerodynamic lift from the SAR antenna — 
amounting to 54 cm. 

Since the original NSWC SEASAT geoid had presumably been in agree- 
ment with the NASA/GSFC SEASAT geoid, the application of the above 
three corrections would lead to a 2 meter discrepancy between NSWC 
and NASA SEASAT geoids (although agreeing with DMAAC). Bill 

Townsend, of NASA Headquarters, asked for a resolution of this 
discrepancy and, in particular, why the Bermuda calibration, per- 
formed with the aid of NASA laser systems, did not detect the 
antenna-feedhorn separation problem. 


ANALYSIS AND RESOLUTION 

For the corrections listed above to the NSWC geoid the following 
assessments have been made: 

1. The Software Problem. There is no reason to question that 
NSWC found a software problem in their analysis program. 
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2. The antenna-feedhorn separation. Indeed the preprocessing of 
the SEASAT altimeter data did neglect the antenna feedhorn 
path length of ~74 cm (MacArthur 1978, Lorell 1979). 

3. The speculation of an aerodynamic lift from the SAR Antenna. 
As has been found by the NASA/GSFC Geodynamics Group and 
the Geodynamics Corporation (Santa Barbara, CA), the calcu- 
lated lift was in error by several orders of magnitude. It 
was later confirmed by NSWC that their calculation was in 
error by the spacecraft mass (a factor of 2000). 

The following questions have been addressed in these analyses: 

1. Should the DMAAC reference geoids be regarded as an accurate 
reference for comparison? 

It has been determined that the 2 meter error found by DMAAC 
was based on a comparison of the NSWC geoid with geoid heights 
of coastal tracking stations obtained from WGS 84. It is 
NASA's assessment that, even if the tracking station heights 
were absolutely correct (and current estimates indicate they 
may be systematically in error by some 40 cm), geoid slopes in 
coastal regions are sufficiently great as to make this test of 
very limited validity. 

2. Why did the Bermuda calibration show the altimeter to have a 
zero bias when the antenna-feedhorn separation had been 
negl ected? 

A thorough review of the Bermuda calibration process 

(Kolenkiewicz and Martin, 1982; Marsh and Martin, 1982; Martin 
and Kolenkiewicz 1981 and Marsh et al., 1982) found no 
problems, either in concept or performance. The calibration 
used data processed the same as that to be operationally 
distributed (with corrections made for altimeter instrument 
delay and for the altimeter offset from the spacecraft center- 
of-mass). The instrument delay applied by JPL was the (APL 
recommended) average of two prelaunch calibrations, three were 
actually made (MacArthur, 1978). Recently, APL has both changed 
their recommended instrument calibration and recomputed the 
antenna-feedhorn offset, with a net change of 58 cm from what 
was applied to the SEASAT data (Kolenkiewicz 1985). As far as 
the inflight calibration is concerned, what was determined was 
the correction to the data as it was actucally preprocessed. 
The conclusion is that the Bermuda calibration did exactly as 
it was purported to do: Calibrate the SEASAT Altimeter . No 

further corrections are necessary to alter the geoid height 
obtained using this calibration, regardless of the number of 
errors that may have been made in processing the altimeter 
data (Kolenkiewicz 1985). 
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FUTURE ANALYSIS 


NASA Personnel are convinced that there is no problem with the 
Bermuda inflight calibration or with SEASAT altimeter data 
processing. This conclusion implies the current state of affairs: 

1. The DMAAC WGS 84 geoid, containing no SEASAT altimetry data is 
2 meters above the NASA/GSFC SEASAT (and GEOS-3) altimetry 
only geoid. 

2. The NSWC SEASAT geoid is 71 to 85 cm above the NASA/GSFC 
SEASAT geoid. 

As requested by NASA Headquarters, offers have been made to assist 
NSVIC in resolving any problems with SEASAT data. Specifically, 
NASA has requested NSWC processed SEASAT altimeter data 

(Kolenkiewicz 1985), which will then be compared with NASA 
processed data. To date, (June 1, 1986), none of the NSWC data 
have been supplied to NASA. 
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CHAPTER III 


OVERVIEW 


GEODESY AND GEOPOTENTIAL 


The development of techniques for determining satellite orbits to 
high precision and accuracy is basic to the extraction of 
geodynamic, oceanographic, gravitational, and other geophysical 
information from observations on satellites. Indeed the develop- 
ment of global gravity models is largely dependent on high preci- 
sion observations to near-earth satellites; the gravity field in 
turn has become a basic data set for understanding the internal 
structure and dynamics of the earth. Similarly, the extraction of 
polar motion and tectonic plate motion from satellite laser ranging 
observations has depended on high accuracy orbit determination. 
Time variations in the orbital elements of Lageos and other satel- 
lites can be related to gravity field and tidal effects, and 
forcing due to interactions of the satellite with electromagnetic 
radiation, charged and neutral particle, and other phenomena. 
Altimetric measurements of the height of a satellite above the 
surface are used for determining the shape of mean sea level (the 
geoid) and the sea surface topography. Such data are then used in 
studies of ocean dynamics. The papers in this chapter focus on 
these elements of geodesy and orbit determi ntaion. The topics dis- 
cussed include the theoretical formulation of the gravita- 
tional perturbation problem for near-ci rcul ar orbits, the Lageos 
orbit, and Lageos orbit perturbations due to various radiative and 
thermal effects. Additionally, several papers are devoted to the 
ongoing geopotential model development for the TOPEX satellite 
(planned for launch in the early 1990's). This oceanographic satel- 
lite will be equipped with a high precision radar altimeter which 
will be used to make precise measurements on the sea surface height. 
The measurements are to be used for studying the dynamics of 
oceanic processes and ocean-atmosphere interactions. To properly 
interpret the data, the TOPEX. orbit must be known with a radial 
error of 10 cm or less. Such a precise orbit can be obtained only 
with an enhanced and improved global gravity field over that which 
is currently available. The relevant papers in this chapter are 
devoted to a discussion of the derivation of this precise gravity 
field through the analysis of high quality satellite observa- 
tions, particularly those involving laser ranging. Finally, an 
analysis of altimetric data over inland sea is presented as an 
example of the utility of radar height determinations to the study 
of tectonophysics processes. 
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The problems in geodesy, geopotential field modeling, and satellite 
orbit determination are strongly coupled to those in Tectonophysi cs 
and Sea Surface Topography. The interested reader is referred to 
the other chapters of this report for further discussions of these 
topics. 

Contributors to this chapter are: Theodore L. Felsentreger, 
Francis J. Lerch, James G. Marsh, David P. Rubincam, Braulio V. 
Sanchez, David E. Smith, and Jean E. Welker 
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LAGEOS ORBIT DETERMINATION 


David E. Smith 


OBJECTIVE 

The orbit of the Lageos satellite has been re-analyzed for the 
period 1983-1985 using improved orbital and earth models, and the 
new Geodyn 2 software system. The analysis was conducted in 
preparation for the next global analysis of Lageos data for the 
determination of tectonic plate motion and crustal deformation. 


BACKGROUND 

Lageos was launched into a high altitude circular orbit in 1976 and 
has been tracked regularly by laser ranging systems ever since. 
Since launch the quality of the tracking has steadily improved from 
the decimeter level to the present day few centimeter level for 
most of the international laser tracking network. Our knowledge of 
and ability to model the many perturbing forces acting on Lageos, 
particularly the earth's gravity field and the ocean tides, has 
also significantly improved over this time. We are now able to 
compute the orbit of Lageos to a much greater precision (and 
probably accuracy) than we could only one or two years ago. Data 
in the period, January 1983 through December 1985, has recently 
been re-analysed in preparation for a new global station coordinate 
solution using our most up-to-date models and software systems. 


RECENT ACCOMPLISHMENTS AND RESULTS 

Thirty-six orbital arcs, each one month in length, have been 
determined and analyzed for the period January 1983 to December 
1985. Each monthly orbital arc contained all the presently avail- 
able Lageos laser data from the international global tracking 
network. The apriori parameter values and models are shown in 
Table 1. The solution to each orbital arc included the orbit, the 
position of the pole and Universal Time (UT), an estimate for the 
value of GM (the product of the earth's mass and gravitational 
constant), and the coordinates of the tracking stations. 

Table 2 summarizes the preliminary solution for each arc. Of 
particular interest in Table 2 are the rms deviations of the 
observations (normal points) from the orbit, which are averaging in 
the 5 to 8 centimeter range, with a few exceptions above 10 
centimeters. These orbital fits are on average several centimeters 
better than were obtained in the SL 6 solution for the same orbital 
arc and reflect the improved models being used in the new SL 7. 
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Because of the adjustments that are made to the station coordinate 
in each arc the rms deviations are indication of overall data 
precision. 

The values obtained for GM are slightly larger in this new solution, 
averaging around 398600.442 km**3sec**-2 for the 36 month period, 
compared to 398600.436 +/- .002 from the 9-year SL 6 solution. At 

the present time it is not known why there has been this small but 

systematic difference from SL 6 but it may be caused by the enlarged 
ocean tidal model used in the new SL 7 solution increasing the energy 
of the system. The new value presently being obtained is more consis- 
tent with the value of GM being obtained from Lageos long orbital arc 
analyses at the University of Texas at Austin, and also from the 

analyses of lunar laser data at NASA/JPL; both of whom are recovering 

values in the 398600.440 to 398600.410 km**3sec**-2 range. Table 3 
shows the values of GM obtained from the three annual solutions for 
the station coordinates. 

Two accelerations of the spacecraft were recovered in each orbital 
arc; each acceleration covering a 15 day period. In previous solu- 
tions only one acceleration had been recovered for the full 30 day 
period. The improvement in the orbital, fit (rms) as a result of 
increasing the number of accelerations from 1 to 2 per arc was very 
small, less than a centimeter in most cases. However, the values of 
the acclerations that were obtained appeared to be well determined, 
providing a historical variation over the three year period very simi- 
lar to that obtained from earlier solutions. In addition, the Lageos 
along-track acceleration accumulates to about 5 meters (on average) 
in along-track position in 15 days and it was felt that this should 
be easily recoverable and would improve the overall accuracy of the 
orbit computation. 


FUTURE EMPHASIS 

The purpose of this analysis has been to prepare for the next (SL7 ) 
solution for global plate motions. The data set described here will, 
over the next few months, be extended back to 1976 and forward to 
include the early months of 1986. At the present time the Lageos 
orbital studies appear to be substantially better in the new analyses 
than earlier ones. In addition, there is clear evidence of the steady 
and continued improvement in the quality of the data throughout the 3 
year period. 
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Table 1. Nominal Parameter Values and Models of SL 7 


Lageos: 

Area 

Mass 

Earth Models: 

Gravity field 

Earth tides 
Ocean tides 

Body tides 
Gravity tide 
Earth rotation 
Planets 
GM 

\ 

Flattening 
Station Coordinates 

Orbit Models: 

Data: 

Arc length: 

Acceleration: 

Solar Radiation Parameter: 

Albedo: 

Relativity: 

Software System: 


0.28274 m 2 
406.965 kg 

GEM L2 (PGS 1715C) 

(CS(2,1) = 0) 

Wahr 

Schwiderski, expanded by Christodoulidis to approx. 600 terms out to 
degree 6 and order 2. 

= 0.609, 1 2 = 0.0852 

1^ = 0.30, kj = 0; phases zero 

Lageos, SL6 polar motion; BIH Universal time 

Gravity fields of Venus, Mars, Jupiter, Saturn 

398600.4359 x !0 9 m 3 sec -2 

6378137.0 meters 

1/298.257 

SL 6, rotated to be consistent with BIH mean pole of 1980-84, epoch 
1983.0, with tectonic motion model AM 1-2 of Minster, Jordan 


2 minute normal points, Herstmonceux definition 
30 or 35 days 

2 per arc (one for each 15 days) 

1 per arc 

None 

None 

Geodyn II E Version 8601.6 
Geodyn II S Version 8601.4 
TDF Version 8601.0 

Solve Version 8604.4 
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Table 2. 



# N. 
POINTS 

RMS 

GM 

SOLRAD 

ACCEL1 

ACCEL2 

JAN 1983 

2156 

.0871 

.436103 

1.1280 

-3.4728 

-1.9100 

FEB 1983 

1854 

.0938 

.436039 

1.1187 

-1.2659 

-1.0149 

MAR 1983 

2360 

.1163 

.447487 

1.1358 

-0.3659 

-0.3357 

APR 1983 

2650 

.1153 

.446485 

1.1581 

-0.3486 

-0.5177 

MAY 1983 

2628 

.1101 

.442030 

1.1747 

-0.6822 

-3.3470 

JUN 1983 

2750 

.0892 

.438976 

1.1754 

-2.9646 

-3.0004 

JUL 1983 

2504 

.0850 

.440372 

1.1710 

-3.1163 

-2.7362 

AUG 1983 

2358 

.0719 

.437685 

1.1615 

-3.3236 

-3.1336 

SEP 1983 

3679 

.0832 

.440812 

1.1490 

-3.3022 

-3.2352 

OCT 1983 

5042 

.0863 

.438748 

1.1342 

-3.5701 

-3.3510 

NOV 1983 

3817 

.0660 

.440942 

1.1360 

-3.9255 

-3.7097 

DEC 1983 

3136 

.0851 

.439168 

1.1344 

-3.8881 

-3.7083 

JAN 1984 

3754 

.0830 

.440461 

1.1381 

-4.2159 

-4.1993 

FEB 1984 

3428 

.0824 

.447869* 

1.1391 

-3.9027 

-3.9029 

MAR 1984 

3030 

.0895 

.442593 

1.1449 

-4.3083 

-3.6693 

APR 1984 

5003 

.1149* 

.447785* 

1.1165* 

-3.6012 

-3.7488 

MAY 1984 

5612 

.0957 

.443547 

1.0884* 

-3.4054 

-3.5003 

JUN 1984 

6007 

.1019* 

.442965 

1.1393 

-3.6195 

-3.5948 

JUL 1984 

6269 

.0991 

.447069* 

1.1357 

-3.6474 

-3.7233 

AUG 1984 

6915 

.0876 

.442681 

1.1290 

-3.2606 

-2.9631 

SEP 1984 

6159 

.0783 

.441776 

1.1298 

-2.4225 

-2.6875 

OCT 1984 

5343 

.0675 

.441947 

1.1274 

-2.2271 

-2.2892 

NOV 1984 

4197 

.0805 

.441734 

1.1269 

-3.8781 

-3.8518 

DEC 1984 

3420 

.0808 

.441107 

1.1338 

-3.9171 

-3.9094 

JAN 1985 

2726 

.0588 

.443669 

1.1274 

-4.0738 

-3.3464 

FEB 1985 

2519 

.0689 

.443165 

1.1633* 

-4.4251 

-3.3776 

MAR 1985 

3125 

.0744 

.440077 

1.1462 

-3.9302 

- 3.7261 

APR 1985 

4012 

.0822 

.443068 

1.1443 

-3.2980 

-2.3004 

MAY 1985 

3223 

.0636 

.440662 

1.1416 

-5.1747 

-6.8523 

JUN 1985 

5074 

.0669 

.439999 

1.1405 

-5.4280 

-5.1037 

JUL 1985 

5542 

.0698 

.440364 

1.1406 

- 4.7332 

-4.7977 

AUG 1985 

5654 

.0627 

.441228 

1.1357 

-4.0265 

-4.0593 

SEP 1985 

7036 

.0642 

.437719 

1.1324 

-3.8560 

-4.2805 

OCT 1985 

5027 

.0608 

.437251 

1.1342 

-3.9852 

- 4.7499 

NOV 1985 

4283 

.0540 

.442727 

1.1364 

-3.3935 

-3.6056 

DEC 1985 

2628 

.0559 

.439182 

1.1363 

-3.4028 

-3.3641 
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Orbital 

RMS 


Table 3. 


GM 


1983 

9.9 cm 

398600.4405 

± 

.0010 km 3 sec -2 

1984 

8.2 cm 

398600.4440 

± 

.0008 km 3 sec~ 2 

1985 

7.1 cm 

398600.4417 

± 

.0008 km 3 sec -2 


95 


GRAVITATIONAL PERTURRATION OF RADIAL POSITION 
AND VELOCITY FOR THE NEAR -CIRCULAR 
SATELLITE ORBIT 


Francis J. Lerch 


OBJECTIVE 

A new method is developed to derive the radial position and velocity 
spectrum of gravitational perturbations for a near-circular orbit. 
These results are compared with those currently obtained from 
linear perturbations of Kepler variables which are near-singular 
for orbits with small eccentricity. 


BACKGROUND 

The radial position and velocity spectrum of gravitational harmonics 
for near-circular orbits are respectively i/seful in the analysis of 
altimeter missions (Tapley and Rosborough, 1985; Wagner, 1985) and 
the NASA planned Geopotential Research Mission (Kaula. 1983: Colombo 
1984; Wagner, 1983). 

Orbits of small eccentricity (e = .003) give rise to non-linear 
Kepler perturbations because e occurs as a divisor in the Kepler 
elements of mean anomaly (M) and argument of perigee ( to) . Kozai 
(1961) and Izsak (1961) derived radial perturbations for leading 
zonal terms of the potential using non-singular variables (e cosu> 
e sinui ). Using these same variables Cook (1966) showed how the non- 
linear effects of long period in e and u> may be stabilized (frozen 
orbit), but short period effects due to J 2 remain. Kaula (1983) 
also suggested to employ these non-singular variables due to the 
small e divisor in GRM orbits. Colombo (1985) used Hill's equations 
to avoid the singularity problem of small e and used the extensive 
zonal field as a reference orbit for GRM analysis. 


ACCOMPLISHMENTS 

In the present analysis non-singular variables 

u = e cos M (1) 

w = e sin M 

are employed which are directly related to the radial distance (r) 
and velocity (r) as follows: 

r = a(l— e cos E) = a(l— e cos M) + Ole 2 ) 

a 2 (2) 

r = — e n sin E = aen sin M + O (e 2 ) . 


I 
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where dominant effects of the perturbations are carried out through 
first order in eccentricity. The variational equations from (1) 
are related to the Kepler form as 

u = e cos M - eM sin M 

(3) 

w = e sin M + eM cos M 

where e and M are given by the Lagrangian equations, as in the 
Theory of Satellite Geodesy (TSG 3.38) by £aula (1966), namely 


. = 1 - e 2 5R _ (1-e 2 ) 5R 
na 2 e 3M na 2 e 5w 


M = n - 


1 - e 2 5R _ 2^3R 
na 2 e 5e na 5a 


where e and M are retained for convenience but are understood by (1) 
to be 

e = (u 2 + w 2 ) 1/2 


M = tan 1 w/u 

The other non-singular Kepler variables are u> + M, a, i, and ft. 

The disturbance potential R is represented in a Fourier series in 
Table 1 with angular arguments 4 1 and amplitudes R. 

Perturbations, Au = u - u 0 and Aw = w-w 0 , use a mean reference orbit 
(TSG 3.74) which is a secularly moving ellijpse with, angular rates 
(u> 0 , ft 0 , n 0 ) and mean elements (a 0 , e 0 , i 0 , ft 0 , uT 0 , M 0 ) as follows: 


Uq C 0 cos M 0 


w„ = e n sin M„ 


M 0 — i^t M 0 , ft 0 o t ft 0 > w 0 co 0 t oj 0 

Perturbative derivatives, Au and aw, are obtained from (3) thru (6) 
to the desired accuracy in e as follows: 

Ail = -n,, Aw + u' 


Aw = ^ Au + w' 


w'l = 

lu'J 


cos M sin M eM' 
-sin M cos M e 


where fi‘ = M - n 0 . 


1 /5R_ _ 3R 

e \5M 5w 
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The eccentricity function of the disturbance R (Table 1) is approxi- 
mated for dominant terms in eccentricity as 


Gfpq(e) = G fpq (0) + e G'fpq(O) 

Gfpo(e) = 1 + 0(e 2 ) for q = 0 

G/ N (e) = | [f + 2q (( - 2p) + 1] + 0(e 3 ) for q = ± 1 . 


The Fourier components of (9) from Table 1 are zero for q=0 

q = + 1 


eM' 

R fmpq 

-Sq 

e 

na 2 

fmpq 

qS'q 


where S q = S/ mpq denote the Fourier terms with arguments 

itmpq = (f-2p) (w+M) + qM + m(O-0) 


and where S' q 



Ffmp(i) G'fpqfO). 


Since 


S q 


[cos *f mpq sin, mpq ] 


Qm 

-Sfm 


f-m 


or 


even 


_ S fm 

- c fm _ 


t — m 
odd 


then 


Sq 


cos qM 

sin qM 

So 

_ S V 


-sin qM 

cos qM 

_So 


Using (11) and (13) with (8) the near-singular variable M may 
minated as follows: 


( 10 ) 

and for 
( 11 ) 

( 12 ) 


(13) 
be el i - 
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cos M sin M 


w 


= A 

fmpq 


-sin M cos M 


-1 O' 
0 q 


cos qM sin qM 
-sin qM cos qM 


J 


S'„ 


= A 



for q - ±1 


( 14 ) 


where A = R' fmpq /na 2 , S G = S, mpo . 

Equation (14) shows that q is eliminated from the arguments in (12) 
where ^ampq becomes I'ampo- 


From equations (11) and (14) the near-singular variables M and e 
are both eliminated from (8) and thus (7) may now be integrated (term 
by term) by evaluating the Fourier components on the reference orbit 
(6) with the mean elements. Here the amplitudes (19) are constant 
and the arguments from (6) and (12) are secular as follows: 

♦O = ^mpo = it + ^ 0 

*o = ifonpo = (f-2p) («o + nj + m (fi 0 - 0) 


*o = hnpo = ( ( - 2p) Go + M„) + m(fl 0 - e o ) . 


The simultaneous system (7) may be integrated for a particular solu- 
tion since the homogeneous solution is given by the reference terms 
(6) for u 0 and w 0 . Perturbational derivatives from (12) are of the 
following form where a and b are constants representing coefficient 
values (C£ m , Sam), namely 


w' 


-A H„ ” 


H 0 


a cos \I/ Q + b sin ^ 0 ~1 

u' 

fmpq 

qA H' 0 


_ H 'o_ 


1 

X ^ 

TJ *0 


and the perturbations have a solution from (7) for q - + 1 as 


Aw 


(k - i^q) H'o 


_ A 


Au 


(qi/'o — n o) Ho 


fmpq 

_ — 1 


( 16 ) 


where \l/ Q is given in (15) above. 
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RESULTS AND SIGNIFICANCE 


Following equation (16) the solutions for the perturbations Au 
and Aw for Ampq components are given in Table 2. From (1) and (2) 
the dominant perturbations in eccentricity for radial position and 
velocity are respectively. 


(17) 

(18) 


A special point of significance is that equation (14) gives 
the same result when derived from linear perturbations for Ae and 
AM (TSG 3.76) and when (14) is evaluated on the reference orbit. 
This is because the near-singular variable M (including non-linear 
effects) vanishes from the equation as a separate variable, at 
least for this case where only the dominant effects in eccentricity 
are being treated. Wagner (1983 and 1985) applied linear perturba- 
tions ( Ae ,AM) to obtain Ar and Ar and his results agree with 
those derived in this report. These results are shown directly 
from the formal solution to (7) and (8) with use of (11) throuqh 
(14) as follows: 

for £mpq components with q = +_ 1: 


Aw cosM 0 sinM 0 

Au -sinM 0 cosM 0 

where on the reference orbit 

w cosM 0 sinM 0 w ' e 0 M' 

u — sinM c cosM 0 u' Ae 

J o L - 1 L J o L J 

hence 

e 0 AM 
Ae 


FUTURE EMPHASIS 

The accuracy of the position and velocity gravitational perturba- 
tions may be improved by extending the theory to include higher 
order terns in eccentricity. 






Alfmpq J Ar f m pq dt 
Afympq an c Aw* mpq 
for q = ± 1 . 
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Table 1. Disturbance Potential R* in Orbital Coordinates With Fourier Components and Related 
Derivatives 



*TSG 3.70 
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THE YARKOVSKY EFFECT, THE SCHACH EFFECT, 
AND LAGEOS 


David P. Rubincam 


OBJECTIVE 

The objective of this study is to find whether the Yarkovsky effect 
or the Schach effect accelerate Lageos at the 10"1 2 ms -2 level. 


BACKGROUND 

The Lageos satellite was launched in 1976 to determine tectonic 
plate motion, polar motion, and length-of-day. Lageos is a spheri- 
cal, completely passive satellite. Its aluminum shell is studded 
with fused silica retroref lectors to permit laser tracking. Inside 
Lageos is a cylindrical beryllium copper core. The orbit is nearly 
circular and inclined to the earth's equator by 110 degrees. The 
altitude is 5900 km. 

Lageos is experiencing a fluctuating along-track acceleration of 
unknown origin of order 3 x 10 -12 ms -2 (see Figure 1). It is 
important to understand the origin of the forces causing the 
acceleration, not only from the standpoint of scientific curiosity, 
but also to aid in the performance of Lageos' mission by modeling 
them. The Yarkovsky and Schach effects have heretofore not been 
examined in detail to assess their importance for Lageos (Rubincam 
1982; Barlier et al . , 1985). Both effects are radiation-reaction 
forces. 

The Yarkovsky effect can be qualitatively understood as follows. 
The satellite spins while direct sunlight impinges on its surface, 
warming it up. Because the satellite has thermal inertia, the 
hottest spot on the satellite is not the subsolar point. There is 
delayed heating, so that the hottest spot is on the "afternoon" 
side of the satellite, just as the hottest time of day on earth is 
the afternoon, instead of noon. 

Lageos emits infrared radiation due to the solar heating, and the 
radiation carries away momentum. There is a net momentum imbalance 
due to the asymmetric temperature distribution, so that Lageos 
feels a force. Moreover, the force is not along the sun-satellite 
line, due to the delayed heating mentioned above. This is the 
Yarkovsky effect (e.g., Opik, 1951; Burns et al . , 1979; Rubincam, 
1982; Barlier et al ., 1985). 


105 



The Schach effect is more subtle; it operates only when the satellite 
enters the earth's shadow. There is an asymmetric north-south tem- 
perature distribution across Lageos' surface due to solar heating, 
which gives a force along the spin axis. The satellite cools off 
slightly as it enters the shadow, then heats up again as it exits. 
This causes the force to vary. Once again there is a thermal lag 
due to the satellite's thermal inertia, causing a net force when 
averaging over one revolution (Rubincam, 1982). 


RECENT DEVELOPMENTS 

The Yarkovsky and Schach effects for Lageos' aluminum have been 
modeled by assuming Lageos is a solid homogeneous sphere of radius 
Rl, density p, thermal conductivity K, and specific heat Cp and 
rotating at angular speed w. The sun was assumed to be infinitely 
far away. The temperature distribution T inside Lageos was found 
by solving the heat conduction equation 

aT 

pCn — = V 2 T 
at 

in spherical harmonics (e.g., Merzbacher, 1970, pp. 194-198) subject 
to the boundary condition 

3T 

eoT 4 + K — = (1-A) F 
ar 

at the surface. Here e is the infrared emissivity, F the irra- 
diance, A the albedo, and a the Stefan-Boltzmann constant. (The 
boundary condition was linearized by assuming T=T 0 + AT, where T 0 
is a constant and T 0 >> AT. The only term contributing to the tem- 
perature asymmetry is the i = 1 term, where i is the degree of 
the spherical harmonic (Rubincam and Weiss, 1986). 

The results for the Yarkovsky effect were as follows. F Qr ]ithe 
reasonable values Ri = 0.3 m, p = 3634 kg nTv, K = 150 Wm" K" 

C = 672 Jkg" 1 K' 1 , E = A = 0.3, and w = 1 rad s' 1 , the amplitude of 
the temperature asymmetry was about 0.01K and the lag angle was 
about 45 degrees (the maximum lag angle). These give an accelera- 
tion of about 0.1 x 10"12 m s“2 f 0 r Lageos. 

For the Schach effect the results were a maximum temperature 
difference of about 0.8K between the poles of Lageos, with a thermal 
lag angle of. about 7.5 degrees. The net along-track acceleration 
is about 0.2 x 10"^ ms - ^. 
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SIGNIFICANCE AND FUTURE EMPHASIS 


The Yarkovsky effect and the Schach effect both give along-track 
accelerations about an order of magnitude below those observed 
(Figure 1). Hence they cannot be the cause of the unmodeled along- 
track force, and the search should continue. 
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PERTURBATIONS ON LAGEOS' ORBIT 
DUE TO INFRARED RADIATION 


David P. Rubincam 


OBJECTIVE 

The objective of this study is to find whether infrared radiation 
pressure from the earth gives long-period accelerations at the 
10' 12 ms* 2 level. 


BACKGROUND 

The rationale for investigating the unknown forces on Lageos are 
given in "The Yarkovsky Effect, the Schach Effect, and Lageos' 
elsewhere in this issue. 

The effect of reflected sunlight on Lageos' orbit when the earth 
reflects like Lambert's law is given by Rubincam and Weiss (1986). 
.They found that perturbations of Lageos ', orbit were negligible, in 
agreement with Barlier et al . (1985). The next logical step was to 
find the effect of the earth's infrared radiation pressure on the 
orbit. The effect was presumed to be even smaller than that due to 
reflected sunlight, due to the way infrared radiation varies 
smoothly across the earth's surface, in contrast to the reflected 
sunlight. 


RECENT ACCOMPLISHMENTS 

The formalism developed by Rubincam and Weiss (1986) has been 
modified to handle infrared radiation emission by the eartiu The 
modifications are: all d £ =0, except dg = 1; F° s (r° s /r s )^ = 1; 

and replace the Anki with E^t, where the Emim are the spherical 
harmonic coefficients of the earth's existence/ so that Eqoi = 232 
Wm , E 2 Q 1 = -26.6 Wm" , etc. These suffice to give the correct 
equations for infrared radiation, assuming the emission follows 
Lambert's law. 

Only the accelerations for the N=0 and N=2, K=0 terms were found 
since these are the two dominant terms in the spherical harmonic 
expansion. The N=0 term provided a check on the modifications, 
since this term gives only a constant radial force. The formalism 
performed as expected, giving the radial force of correct magnitude 
and zero for the along-track and normal -to-orbit forces. 

The N=2 or E 20 L term 9 ave only short-period radial accelerations of 
order 14 x 10"* 2 m s~ 2 ,and along-track accelerations of order 6 x 
lO - * 2 m s _2 . There were no long-period perturbations in the semi- 
major axis to order e in the eccentricity. 
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SIGNIFICANCE AND FUTURE EMPHASIS 

Infrared radiation pressure from the earth cannot explain the 
fluctuations in along-track acceleration of order 3 x 10" 1Z m s"^ 
observed on Lageos. This was expected and agrees with Sehnal 
(1981). Hence the search for the origin of the unmodeled along- 
track forces must turn elsewhere. 
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GEOPOTENTIAL MODEL DEVELOPMENT FOR TOPEX 


James G. Marsh 


OBJECTIVE 

The TOPEX oceanographic satellite mission which is planned for the 
early 1990 ' s will carry as a prime instrument a radar altimeter 
system. In order to use altimetric data for the mapping of the ocean 
surface, radial orbit information with an accuracy approaching a 
decimeter is required. To achieve this level of accuracy, at least 
a factor of 2 improvement in the modeling of the earth's gravity 
field is required. A major effort is underway in the Geodynamics 
Branch for the computation of a new earth gravity field model which 
will support the TOPEX precision orbit determination requirements. 


BACKGROUND 

Satellite altimeters measure the height of the satellite above the 
sea surface. When these measurements are combined with the satel- 
lite position information determined from the tracking observations 
of the satellite, the sea surface height with respect to a refe- 
rence earth ellipsoid can be computed. The information on temporal 
and spatial variations in sea surface height will lead to a substan- 
tial improvement in our understanding of the ocean circulation and 
global ocean processes. The satellite altimeter hardware system 
planned for TOPEX will have an accuracy of about 2 cm and the total 
system error budget (including the radial orbit error) is on the 
order of 10 cm. In order to meet this stringent requirement, a 
significant improvement in the model of the earth's gravity field 
is required. In support of this research a comprehensive set of 
computer programs has been developed. Most of this software was 
developed for use on the IBM 360 series computers. The recent 
availability of the Cyber 205 Vector Processing Super Computer has 
provided a significant enhancement in the capability to perform 
such massive numerical computations. For example, in the actual 
least squares solution for the gravity coefficients, a factor of 30 
improvement in speed has been achieved. 

Significant improvements in the quality of laser tracking data have 
been achieved since the computation of the previous general GSFC 
gravity model GEM 10B. These new data form a major portion of the 
observational data base used for the new solution. 
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RECENT ACCOMPLISHMENTS 


The basic objective of the TOPEX Gravity model activity within the 
Geodynamics Branch during the past year has been the computation of 
a preliminary gravity model solution based upon a set of the most 
important satellite tracking data. 

The total observational data base available for analysis consists 
of observations on 60 satellites with a range of accuracies from 
tens-of-meters to a few cm. In order to establish an optimum data 
base for observations on the preliminary solution, a comprehensive 
satellite ranking system was established which considered such 
factors as unique orbit contribution, data accuracy, and non- 
conservative force model errors. A total of 11 satellites was 
selected. The recently acquired laser data on the geodetic space- 
craft formed an important part of this data base. Orbital arcs, 
generally 5 to 7 days in length (30 days for Lageos) were selected. 
An extensive data validation effort resulted in a set of relatively 
error-free data arcs for each satellite. In this preliminary 
solution the coordinates of the tracking stations were derived 
using datum shifts based upon the best Lageos solution and were 
held fixed. Solid earth and ocean tides produce significant pertur- 
bations on the satellite orbits considered for the analyses. An 
extensive long wavelength earth and ocean tidal model has been 
implemented for the data reduction. In order to insure uniformity 
in the orbit computations for all 11 satellites a standard set of 
constants has been adopted. 

The total software system for orbit and geodetic parameter estima- 
tion has been converted to the Cyber 205 Vector Processing Computer. 
Comprehensive data analysis and matrix generation procedures 
compatible with the Cyber resources have been created. The availa- 
bility of the Cyber 205 computer system has permitted a far 
superior data analysis capability than was possible on the earlier 
IBM 360 system. 

The analyses to date have resulted in the computation of a prelimi- 
nary gravity model based upon over 500 arcs of data on the 11 
satellites. The model is complete to degree and order 36 and 
is based upon the best available set of constants and force 
models. In addition, a version of the model has been produced 
which has permitted the adjustment of some 66 ocean tidal terms. 

Accuracy checks of the model based upon rms fits to tracking data 
have demonstrated significant improvement for all satellites even 
when the results are compared with earlier "tailored" models. For 
example, on Starlette the rms fits were reduced from over 70 cm for 
an earlier Starlette tailored model (PGS-1331) to 30 cm for the new 
model . 
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SIGNIFICANCE 


In order to meet the radial orbit accuracy requirement of the TOPEX 
Project, a factor of at least 2 improvement in the gravity field is 
required. The preliminary results computed to date are indicative 
that substantial progress can be made in the reduction of the 
gravity field errors. The overall gravity field improvement effort 
will be carried out over a time period of several years with the 
final model being available prior to the launch of the TOPEX 
satel 1 ite. 


FUTURE EMPHASIS 

Future emphasis will be directed toward the incorporation of 
additional data sets such as, altimetry, satel 1 ite-to-satel 1 ite 
tracking and a re-assessment of the value of the lower priority data 
sets. Surface gravity data will be incorporated in order to provide 
better definition of the higher degree and order coefficients. A 
calibrated set of statistics will be derived which accurately 
describe the errors in the new gravity field so that the performance 
for TOPEX orbit computations can be assessed. 
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ANALYSIS OF GEOS-1 LASER RANGING DATA FOR TOPEX 


Theodore L. Felsentreger 


OBJECTIVE 

The intent of this work is to prepare a 2-year set of GEOS-1 laser 
ranging data for inclusion in a data base which is to be used in a 
solution for the gravity field of the Earth. This new gravity 
field model is needed so that the orbit of the proposed TOPEX 
satellite can be determined accurately enough for the mission to be 
accomplished. 


BACKGROUND 

The TOPEX (Ocean TOPography Experiment) is a dedicated satellite 
whose primary instrument will be an altimeter with an anticipated 
precision of 2 cm. This altimeter will be used to measure sea 
surface elevation to an expected accuracy of a few centimeters. 
This accuracy is necessary for worthwile information about ocean 
circulation to be obtained. Thus, the accuracy of the TOPEX orbit 
is critical, and this implies the use of an accurate gravity field 
model for the Earth in the orbit determination procedure. 

A large data base of satellite tracking data exists at GSFC. An 
effort is underway to collect and process this data (consisting of 
Optical, Doppler and Laser) for use in the derivation of a model 
for the Earth's gravity field. Laser tracking data for the GEOS-1 
satellite is part of this data base. 


RECENT ACCOMPLISHMENTS 

GEOS-1 laser data from the period January 20, 1977 to December 14, 
1978 have been chosen for analysis. This time span more than covers 
one period of the argument of perigee, thus providing good global 
coverage. The data involves both SAO and NASA stations. 

The first step in the procedure was to catalog the data and divide 
it into 5-day arcs, eliminating those time periods with little or no 
coverage. Attention was given to the number of passes and the 
number of stations involved in any 5-day period. A total of 104 
arcs survived this scrutiny. Table 1 provides a summary of the 
satellite's orbit and the tracking data. 


Receding page plank not filmed 
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The NASA data was at a frequency of 1 measurement/sec., with 1 
measurement/2.5 sec. for the SAO data. It was decided to choose 
every third NASA observation and every SAO observation to provide a 

more even weighting of the data. Then, using estimates of the 

position and velocity vectors of the satellite for the arc epochs, 
nominal values for air drag, solar radiation pressure and solid- 
earth tidal parameters, and the GEM-10B gravity field, the arcs 

were converged. In the convergence process, the position and velo- 
city vectors, air drag and solar radiation pressure parameters were 
adjusted for each arc. The purpose of the convergence is twofold: 
(1) to obtain more accurate position and velocity vectors prepara- 
tory to the creation of the matrix of normal equtions ("E"-matrix) 
to be used in the gravity field solution, and (2) to identify and 
delete nonreliable measurements and/or passes. One air drag coeffi- 
cient (Cq) for each day of a 5-day arc and one solar radiation 

pressure (Cr) coefficient for the whole arc were solved for. The 
data reduction was performed using the GEODYN computer program. A 
total of 53 arcs have been processed, covering the time period 
2/1/78 - 12/14/78. The results are given in Table 2. An rms value 
for each 5-day arc has been computed — this provides an indication of 
the overall fit to the data for the arc. The rms values are gene- 
rally in the 1-2 meter range — this is quite good, considering that 
there is a resonant effect on the orbit caused by the 12th order 
coefficients in the Earth's gravity field. 


FUTURE EMPHASIS 

Future activities will concentrate on processing the data for the 
rest of the arcs. Also, the effects of the ocean tides will be 
modeled and a more extensive solid-earth tide model will be used. 
Closer attention will be paid to eliminating nonreliable measure- 
ments from the data. Finally, the matrices of normal equations 

will be formed to be included in the solution for the gravity 
fi eld. 
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Table 1. Data Summary 


i 


• SATELLITE : GEOS-1 

• ORBITAL ELEMENTS : 

Semimajor axis = 8080 km Inclination = 59°.4 
Eccentricity = .07 Period of arg. perigee ~ 540 days 

• TIME PERIOD : 1/20/77 - 12/14/78 

• ARC LENGTH : 5 Days 

• DATA : SAO + NASA LASER 

• NO. ARCS (INCL. NASA) : 104 (60) 

• NO. OBS. : ~ 155,000 

• NO. PASSES (MIN-MAX): 2274 (13-33) 


I 
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Table 2. GEOS-1 Orbit Determination Summary 



Arc Epoch 
YYMMDD 

No. Obs. 

RMS (m) 


780201 

1091 

2.9030 

* 

209 

1052 

1.3375 

* 

217 

1290 

2.4337 

* 

222 

1666 

1.7366 

* 

308 

876 

1.1719 

* 

314 

1031 

1.3511 

* 

322 

836 

1.1865 

* 

330 

935 

1.1273 

* 

404 

974 

1.5250 


413 

892 

1.1264 


419 

954 

1.3172 


424 

1512 

1.2363 


429 

977 

1.5580 


504 

1331 

2.0774 


509 

1832 

1.1245 

* 

514 

1055 

1.5260 

* 

520 

1085 

1.5142 


528 

1099 

1.2058 


602 

1468 

1.2645 


607 

1698 

1.1120 

♦ 

613 

1565 

1.0501 

* 

620 

1036 

1.0907 

* 

625 

1493 

1.2779 

* 

630 

1503 

1.2411 

* 

705 

1691 

1.6425 

* 

710 

1563 

1.2306 

* 

715 

1224 

1.2703 

* 

720 

1036 

1.4015 

* 

725 

962 

1.6927 


730 

1694 

1.5088 


804 

1397 

1.5396 

* 

809 

1016 

1.8366 

* 

820 

901 

1.7908 

* 

825 

1012 

1.3705 

* 

830 

1000 

1.1439 

* 

906 

814 

0.9754 

* 

919 

1802 

1.2879 

* 

924 

1339 

1.0878 

* 

929 

1474 

0.8966 

* 

1004 

1666 

2.2681 

* 

1009 

2041 

1.4950 

* 

1014 

1911 

1.5439 

* 

1019 

1211 

1.7310 

* 

1024 

2074 

1.7677 

* 

1029 

1280 

1.5958 

* 

1105 

1239 

1.3640 

* 

1110 

1241 

0.7669 


1115 

1391 

1.0669 

* 

1120 

1588 

1.1684 

* 

1125 

1009 

0.9963 

* 

1204 

1378 

1.0902 

* 

1209 

920 

0.9503 


c D 

Cr 

Arg. Per. Range 
(deg.) 

2.373-4.188 

1.020 

189.3 193.3 

2.296-2.959 

1.269 

194.8-198.9 

0.852-2.203 

1.690 

200.7—203.5 

2.320-3.290 

0.780 

203.5-207.2 

2.299-2.613 

0.850 

212.9—215.9 

2.079-2.598 

1.295 

216.5—220.6 

2.242-2.451 

1.259 

222.4-225.3 

2.311-2.601 

1.198 

227.7—230.5 

1.664-2.328 

1.686 

230.5—234.0 

2.286-2.663 

1.009 

236.7—239.9 

2.391-3.144 

1.360 

240.8—244.1 

1.257-2.133 

1.089 

244.1—247.2 

2.256-3.247 

1.633 

247.2—251.0 

2.474-4.522 

1.586 

251.0-253.5 

2.099-2.613 

1.315 

253.5-257.6 

2.305-2.943 

1.628 

257.6-260.2 

2.346-2.941 

1.628 

261.0-264.9 

2.251-2.668 

1.394 

266.5-269.9 

2.381-2.965 

1.248 

269.9-273.4 

1.805-2.258 

1.316 

273.4—276.5 

1.42^—2.241 

1.377 

277.4-280.3 

1.795-2.302 

1.712 

281.5-285.4 

1.780-3.288 

1.512 

285.4-288.3 

2.332-3.374 

1.284 

288.3-291.8 

1.114-2.339 

1.397 

291.8-295.2 

2.117-2.720 

1.277 

295.2-298.6 

2.282-2.887 

1.138 

298.6-301.8 

1.997-2.293 

1.396 

301.8-304.9 

1.977-2.298 

1.610 

304.9-308.5 

2.261-2.586 

1.114 

308.5-311.2 

2.023-2.330 

1.453 

311.2-315.4 

2.038-2.313 

1.650 

315.4-318.1 

2.304-2.488 

0.805 

322.7-325.3 

2.243-2.373 

1.287 

315.3—329.0 

1.931-2.284 

1.258 

329.0-331.7 

2.003-2.306 

0.963 

333.0-337.1 

2.461-3.545 

1.291 

342.5-345.3 

2.236-2.899 

1.279 

345.3-348.6 

2.239-2.312 

1.269 

348.6-351.9 

2.406-4.418 

1.508 

351.9-354.6 

2.517-4.867 

1.233 

354.6-358.7 

2.840-5.302 

1.342 

358.7- 1.6 

2.629-4.126 

1.696 

1.6- 5.4 

2.434-4.895 

1.829 

5.4- 8.0 

2.387-3.889 

1.676 

8.0- 11.7 

2.330-3.365 

1.481 

12.9- 15.4 

2.130-2.632 

1.393 

15.4- 19.5 

1.904-2.482 

1.257 

19.5- 22.3 

2.272—3.477 

1.166 

22.3- 26.3 

2.300-3.401 

1.253 

26.3- 28.9 

2.404-3.428 

0.547 

31.7— 34.8 

2.317-2.644 

0.929 

34.8- 37.7 


*Arc includes NASA data 
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GEOS-3 ORBIT DETERMINATION 


Braulio V. Sanchez 


OBJECTIVE 

The objective of this work is to determine the epoch state vectors 
for 5-day arcs of laser range data for the GEOS-3 earth orbiting 
satellite. These arcs will be used as part of the data base 
leading to the derivation of a gravity field model for the TOPEX 
Project . 


BACKGROUND 

The Geodynamics Earth and Ocean Satellite, GEOS-3, was launched on 
April 9, 1975. The satellite characteristics and the nominal 

orbital parameters are the following: 

Area: 

Mass: 

Eccentricity: 

Incl i nation: 

Perigee Height: 

Apogee Height: 

Orbital Period: 

Argument of Perigee Period: 

The available data were obtained by both NASA and SAO laser tracking 
stations during the years 1975 and 1976. It is distributed as 
fol lows : 

1975: 196916 meas. 

1976: 193405 meas. 

Total: 389421 meas. (SAO: 18%). 

Past experience indicates that a 5 to 7 day arc length is optimum 
for geodetic satellites at 800 to 1000 km orbit heights. This time 
span provides strong dynamic content without incurring excessive 
nonconservative force effects such as atmospheric drag and solar 
radiation pressure. A 5-day arc for GE0S-3 covers approximately 
the period of the effect produced by the resonant 14th order 
coefficients of the Earth's gravitational field. This effect can 
reach magnitudes of 150 meters in the along-track component. The 
gravitational field used in the computations is the GEM-10B model 
full to degree and order 36, derived from satellite tracking data, 
surface gravity and altimetry. The atmospheric density is obtained 
from the Jacchia 1971 model. 


1.4365 m2 
345.909 kg 
0.0014 
115° 

840 km 
860 km 
102 minutes 
1039 days 
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The data reduction is performed by means of the GEODYN least squares 
parameter estimation and orbit determination system. The GEODYN 
system has the capability to estimate the state vector (position 
and velocity of the satellite) as well as many other kinematic and 
dynamic parameters including parameters appearing in the force 
models such a atmospheric drag coefficients (Cq) and solar radiation 
pressure coefficients (Cr). 


RECENT ACCOMPLISHMENTS 

A total of 41 arcs have been processed, covering the time period 
from May to December of 1975. At this point the only editing 
applied to the data has been the dynamic editing inherent in the 
iterative process of GEODYN. Besides the position and velocity of 
the satellite, epoch values for Cp and Cr have been estimated also. 
The trajectory generated using these estimated parameters is used 
to compute an rms value for each 5-day arc, which provides an 
indication of the overall fit to the data for each arc. The results 
are given in Table 1 below. 


FUTURE EMPHASIS 

Future activities will include the processing of the data available 
for 1976 as well as the incorporation of more sophisticated force 
models for the atmospheric drag. The data will be scrutinized more 
closely in order to eliminate nonreliable measurements. Also the 
effects of ocean tides will be modeled and a more extended earth tide 
model will be used; in generating the results up to now no ocean 
tide effects have been included and the earth tide effects have 
been modeled through k 2 and only. 

The converged arcs will be used to compute the coefficients of the 
matrix appearing in the normal equations. Once the new gravity 
field has been computed, it will be used to process the data again 
and the results will be studied. 
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Table 1 . GEOS-3 Orbit Determination Results 


ARC EPOCH 

NO. OF MEAS 

750511 

2312 

750519 

1087 

750524 

1119 

750603 

1054 

750609 

1179 

750614 

3179 

750619 

2940 

750624 

2273 

750629 

3834 

750704 

3126 

750709 

4449 

750714 

2956 

750719 

4173 

750724 

3317 

750729 

3279 

750803 

4049 

750808 

2942 

750813 

2607 

750818 

3502 

750823 

4484 

750828 

5579 

750902 

3967 

750907 

4835 

750912 

1818 

750919 

1813 

750924 

1733 

750929 

1652 

751004 

1646 

751009 

1630 

751014 

1387 

751019 

1723 

751029 

2348 

751103 

2017 

751108 

1868 

751113 

1166 

751118 

2778 

751123 

2698 

751202 

2299 

751209 

2464 

751216 

6322 

751227 

1096 


RMS 

(meters) 


Cr 

1.425 

3.32 

2.11 

0.663 

4.01 

2.52 

0.504 

3.51 

2.41 

1.834 

3.41 

1.37 

1.462 

2.59 

1.59 

1.493 

2.65 

1.62 

1.551 

2.79 

1.30 

1.510 

2.92 

1.46 

1.670 

3.01 

1.55 

2.160 

3.11 

1.37 

2.175 

3.34 

1.76 

1.301 

2.67 

1.71 

1.745 

2.78 

1.72 

1.557 

2.91 

1.25 

1.019 

3.08 

1.52 

1.548 

2.66 

1.66 

1.222 

2.90 

1 .50 

1.283 

3.34 

1.52 

1.591 

3.16 

1.93 

1.637 

3.83 

1.65 

1.141 

3.83 

1.65 

1.001 

4.62 

1.63 

1.195 

4.00 

1.45 

1.004 

3.31 

1.42 

2.295 

2.21 

1.42 

3.006 

2.34 

0.96 

2.105 

2.73 

1.30 

2.624 

2.50 

1.06 

2.639 

2.42 

1.28 

2.953 

2.79 

1.43 

1.657 

2.95 

1.33 

3.033 

2.54 

1.52 

0.618 

2.30 

1.19 

0.248 

2.77 

1.17 

0.161 

2.69 

1.16 

2.688 

2.67 

1.75 

0.587 

2.28 

0.77 

1.955 

2.54 

1.48 

2.696 

3.72 

0.31 

2.329 

3.06 

1.95 

1.449 

3.44 

0.92 
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A SURFACE CONTOUR OF THE CASPIAN SEA FROM SATELLITE ALTIMETRY 


Jean E. Welker 

I 

■ OBJECTIVES 

There are two aspects of this work: 1) to combine a surface contour 
mapping of the Black and Aral Seas with the current Caspian Sea 
mapping for a regional tectonic analysis; 2) to modify and program 
| some existing mathematical models in order to interpret the alti- 

| metric mapping data over these seas. 


BACKGROUND 

The northern and southern portions of the Caspian Sea have differen- 
ces in their origins, and can be divided into two distinctly 
different tectonic zones. The boundary of this division has been 
defined by the shallow east-west zone which cuts across the Caspian 
Sea at the Apsheron Peninsula, forming -the boundary between the 
northern and southern Caspian depressions. This shallow zone or 
boundary also lies along the projected path of the spine of the 
greater Caucasus, were it to continue and pass directly under the 
Caspian Sea. The formation of the northern Caspian depression has 
been linked to the Paleotethyan ocean closing. The southern 
Caspian depression has been related to an enlarged ocean basin 
resulting from motion along the fault system formed by the same 
Paleotethyan ocean closing (Sengor, 1984; Dercourt, et al., 1986). 
From a broader perspective, the region of the greater Cacasus 
(located between the Black and Caspian Seas) is part of the Alpine- 
Himalayan system of orogenic belts which formed from the double 
closing of the Paleotethys and Neo-tethys oceans beginning in the 
late Triassic. This region has a complex of active tectonic fea- 
tures as shown in Figure 1. In both the northern and southern 
portions of the Caspian Sea, linear tectonic features are either 
projected crossing the sea, or the same feature is shown on either 
side of the sea. The continuation of these features in easterly 
and westerly directions intersect the Black and Aral Seas, forming a 
complex of tectonic features which run into or through the three 
seas. 

In order to shed some light on the tectonic features described 
above, a modified approach has been developed for processing alti- 
metry records over inland seas to obtain surface contours for those 
sea surfaces (Zwally, 1983; Martin, 1983; Brenner, 1983). 
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Altimetry records from both the SEASAT and GEOS-3 satellites have 
been processed for all three seas in order to obtain surface 
contours. With these contours, inferences can be made about the 
type and shape of the tectonic structures beneath the seas. 

In a separate but related effort, an existing two-dimensional 
gravity modeling program is being expanded to three dimensions. 
Inherent limitations in the application of certain boundary 
conditions, such as the length and depth dimensions of tectonic 
structures, are being investigated. The program in its final form 
will be used to model the tectonic configurations beneath the seas 
which have produced the anomalies in the sea surface contours. 


RECENT ACCOMPLISHMENTS 

In order to meet the mapping requirements, an investigation of the 
surface topography of the Caspian Sea was initiated. The satellite 
altimetry data from GEOS-3 and SEASAT provided a data base for the 
measurement of the surface height and structure of the Caspian. 
The concept of calculating surface contours from altimetry data had 
al ready yi el ded results in the open ocean areas (Marsh et al., 
1980, Brown, et al , 1983), and over the continental ice sheets 

(Zwally et al., 1983). Because the reduction of the altimetry data 
over the Caspian Sea required many similar techniques to those 
developed for over the oceans and continental ice sheets, a contour 
mapping of the Caspian was possible. The basic technical obstacle 
to the contour mapping of an inland sea such as the Caspian was the 
smooth gridding and contouring of the unevenly distributed alti- 
metry data, which was especially uneven near the water-land inter- 
face. A similar problem had occurred with the unevenly distributed 
altimetry data over the continental ice sheets. When the tech- 
niques developed for the ice sheets were applied to the Caspian Sea 
data, results were obtained near the shore without degrading the 
results in the central regions of this inland sea. The ground 
tracks from the SEASAT and GEOS-3 satellites and the contour map of 
the Caspian are shown in Figures 2 and 3, respectively. The sur- 
face contour of the Black Sea has also been completed, but has not 
been processed in its final form. The surface contouring of the 
Black Sea has been straight forward, and no major problems have 
been encountered. The number of crossover points from intersecting 
ascending and descending satellite groundtracks has been sufficient 
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to satisfy the surface contouring requirements. An initial paper 
describing the results of the contouring of the Caspian and Black 
Seas and interpretation is in preparation. 

The Aral Sea has presented a problem because of the limited number 
of crossover points from the intersection of ascending and des- 
cending satellite groundtracks. Its general shape is roughly 
circular, with a diameter of approximately 250 km. Crossover 
intersections are scarce and some small islands in the sea have 
contributed to a further reduction in the crossovers. If a proper 
contour of the Aral Sea cannot be constructed, surface profiles for 
the individual satellite tracks will be drawn. 


SIGNIFICANCE 

The successful contour mapping of the Caspian and Black Seas to 
date has proven the viability of a new technique for the investiga- 
tion of inland seas. The mapping of surface contours has provided 
insight in the relation of tectonic features bordering the Caspian 
and Black Seas to one another. With the successful development of 
models to utilize the surface contour data, additional insight into 
the underlying tectonic structure, and the origin and formation of 
these seas may be possible. 

FUTURE EMPHASIS 

A study of the region of the greater Caucasus will continue with 
the development of the model technique discussed above. Lake 
Baikal, in Central Soviet Siberia, is also under study. An assess- 
ment of other inland seas which can be contoured using SEASAT and 
GEOS-3 altimetry is continuing. These include the eastern 
Mediterranean, the Lake Victoria region in southeast Africa, and 
the Great Lakes region on the U.S.-Canadi an border. 
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CHAPTER IV 


PLANETOLOGY 


OVERVIEW 


Although individual members of the Geodynamics Branch have pre- 
viously been interested in problems in planetology, 1985 marked the 
beginning of a concerted effort to establish a major planetary 
research program within the branch. Much of the effort is directed 
toward geodetic studies of the planet Mars through gravity field 
mapping and altimetric studies. The work is a direct extension of 
the well-established expertise and experience within the branch in 
developing gravity field models, performing orbit determinations, 
and analyzing altimetric data for earth satellites. The Mars 
Observer Mission, planned for the early 1990's gives new emphasis 
to these studies. The interest of branch personnel in planetology 
has been further heightened by the applicability of their geophysi- 
cal models of lithospheric deformation and solid-earth dynamics to 
processes on other terrestrial-like planets. 

Contributors to this chapter are: B. Fong Chao, Francis J. Lerch, 
David P. Rubincam and Maria T. Zuber. 
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ANALYSIS OF MARS GRAVITY FIELD FROM 
RADIO TRACKING OF MARS OBSERVER 


Francis J. Lerch 


OBJECTIVES 

The objective of this work is to estimate the accuracy and resolu- 
tion of the gravity spectrum from the Mars Observer (MO) mission 
with earthbased DSN Doppler tracking. A global potential of spheri- 
cal harmonics is modeled with an a priori power spectrum based upon 
earlier DSN tracking of the Viking and Mariner missions. Gravity 
anomaly wavelength resolution is derived from the detailed MO 
accuracy spectrum of harmonics. This information is useful for 
understanding the physics (stress state, density contrasts, flexural 
response, lithospheric thickness, etc.) of the interior of Mars and 
its surface. In conjunction with existing topography and improved 
topography from the MO altimeter it will serve to determine bounds 
for theories and models of isostatic compensation of the Martian 
crust. Also the accuracy of the low degree field will indicate if 
temporal changes can be detected that correspond to seasonal 
variations in the Martian polar caps' and atmosphere. 


BACKGROUND 

The representation for the Martian (Ares) gravitational potential 
at the altitude of the MO spacecraft can be defined in spherical 
harmonics as follows: 


V a (r) 

where 


N i 

~‘T I flu 

r 1‘0 


- GM 


t f I*\ l 

n=0 \ r) 


P £m [C £m cosmX+S l(TJ sinmA] (!) 


GM a 

r a 

p £m 


N 


- position vector of the satellite with coordinates 
(r, 4 >,X,) Areocentric distance, latitude, and 
longitude 

- Areoidal central gravitational constant 

- mean radius of the reference ellipsoid of Mars 

- normalized associated Legendre functions of 
degree i and order m 

- normalized spherical harmonic coefficients which 
are to be estimated to define the gravitational 
model 

- maximum degree representing the size (resolution) 
of the field. 
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N is expected to be as high as 50, equivalent to a horizontal 
resolution of about 200 km. 

Harmonic coefficients have been estimated (Balmino et al., 1982) 
from the previous Martian DSN data complete through only degree and 
order 18 (N=18). These coefficients have been shown (Balmino et 
al., 1982) to have an rms size per degree of 

o(C,S) £m = 13 x 10‘ 5 /£2 (2) 


which is 13 times the size of the coefficients for the Earth as 

estimated by Kaula (1966). Balmino's rule (2) for the coefficient 
size will be used when estimating the harmonic perturbations in 
satellite position and velocity of MO. 

The Mars Orbiter will be launched (1990 planned ) in a near- 

circular orbit at 350 km altitude with a mission lifetime of 1 

Martian year (687 earth days). Tracking data is scheduled from the 
earth-based Deep Space Network of NASA with accuracy of ±0.03 

cm/sec over a 1 minute averaging interval as in the previous 
Mariner 9 mission (Balmino et al., 1982). Using this figure as a 
guide and Balmino's rule (2) for Mars, a spectrum of significant 
perturbations in the velocity of MO orbit can be seen in Figure 1 
from the harmonics through degree 50. 

Velocity component perturbations (radial, crosstrack, and along- 
track) of M0 are summarized in Figure 2 for the same spectrum of 
harmonics. The M0 Doppler signal of gravity would be composed of 
the projection of these components on the line of sight from the 
tracking station. Separability of the harmonics in the solution is 
discussed below. From Figure 1 the sensitivity to a single harmonic 
perturbation on M0 is on averge about 0.01 cm/sec at degree 50, but 
for the root-sum-of-squares of all orders there is approximately 0.1 
cm/sec, considerably higher than the noise level of the data. 


ACCOMPLISHMENTS AND RESULTS 

An error spectrum (Figure 3) was derived for an M0 gravitational 
model of spherical harmonics complete through degree and order 50. 
The spectrum shows an order of magnitude improvement over current 
knowledge (Balmino et al., 1982). The analysis assumed DSN Doppler 
data of 0.3 cm/sec noise for 1 minute data intervals oyer 70 days 
of MO tracking. For simplicity, an average component (Figure 3) of 
alongtrack, radial, and crosstrack velocity was used to represent 
the DSN Doppler tracking coverage (described in Figure 5 of the 
Project Overview of Mars Observer Proposal Information Package, 
0SSA2-85, Volume IV). From the same reference the analysis employed 
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the "frozen orbit" geometry (o> fixed at -90°) of MO along with the 
repeating groundtrack of 56-day period of 30 km uniform spacing. 
These conditions remove the long-term perturbation effects and 
permit the gravitational signal to repeat itself. As a consequence, 
the cross terms of the normal matrix of the velocity components of 
the gravitational signal (Colombo, 1981; Wagner, 1983) nearly 
average to zero (barring small nonlinear effects) except for the 
terms of order m and the same degree parity. This allowed the 
matrix to be simplified to a sparse form of block diagonal terms 
corresponding to each order m. Employing the Fourier representation 
(Kaula, 1966) of the potential perturbations, the signal was 
averaged analytically over the span of the repeat cycle, thus 
simplifying point-by-point computations. Because of the occulta- 
tion of the radial and velocity components as given in the Project 
Overview (Figure 13, Volume IV), crossterms will not completely 
vanish. 

From Figure 3 of the error spectrum, the cross-track effects are 
more dominant in the solution than radial or along-track effects. 
Cross-track effects are not occulted by Mars when the orbital plane 
is normal to the line-of-sight. About 1.9 years of tracking cover- 
age (687 days) are available and about 1/10 the quantity of data 
(70 days) is assumed for the solution. By calibrating the error 
spectrum from actual but spotty satellite tracking coverage from 
Earth satellites, the MO error spectrum (Figure 3) was degraded by 
a factor 6 from the ideal solution. Combining this degradation 
with that from the data (/10), a total reduction of a factor of 
19 from the ideal value has been employed as a conservative esti- 
mate of the accuracy for the MO model. The solution also assumed 
use of £ pri ori knowledge (Figure 3) for the gravitational field 
from previous Martian missions. If this a priori assumption is 
removed, there is essentially no change in the error spectrum, 
except at high degrees where the error curve departs along a tangen- 
tial line at degree 30 and approaches closely to the a priori 
estimate at degree 50. 


SIGNIFICANCE AND FUTURE EMPHASIS 

The error spectrum of the harmonics was used to estimate the global 
rms mean gravity anomaly errors for various block sizes (half 
wavelength) as indicated in the table below. The results for the 
400 km blocksize, showing a 6.9% error of total power, will be 
needed to estimate useful density contrasts in the Martian crust 
(e.g., basalt/granite). Wavelengths shorter than 1000 km (500 km 
half wavelength) are expected to show negligible flexural response 
(no compensation). Therefore, at these wavelengths, the gravity 
anomalies should result from the near-surface density contrasts in 
the crust. 
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MARS MEAN GRAVITY ANOMALY* 


POWER SPECTRUM - 13X10' 5 /*. 2 COEFF. MODEL 
MO MODEL - 50X50, 70 DAYS DSN DATA, 0.03 cm/sec. 


MEAN ANOMALY 
BLOCK SIZE 
(1/2 Wavelength) 

TOTAL 

POWER* 

M0 MODEL 
COMMISSION 
ERROR* 

M0 MODEL 
TOTAL ERROR* 

TOTAL PERCENT 
ERROR (%) 
Error/Power 

1.0° 


132.1 mgals 

11.0 mgals 

72.4 mgals 

54.8 

3.5° 

(200 km) 

107.7 

8.5 

26.4 

24.5 

5.0° 

(300 km) 

99.8 

6.0 

13.0 

13.0 

O 

O 

• 

(400 km) 

91.9 

2.8 

6.3 

6.9 

10.0° 


82.9 

1.3 

4.7 

5.7 

♦Gravity 

anomaly 

averaged over a 

block. 




Density structure may be analyzed from the gravity model in 
conjunction with the topography obtained from the MO altimeter and 
will increase our understanding of the Martian lithosphere and its 
evolution. 

Modeling of the MO orbit and gravitational signal from earth-based 
DSN tracking will be developed in order to recover the gravity 
field to the accuracy specified in Figure 3. This developmental 
program will be tested with existing DSN data from Viking and 
Mariner missions. 
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HARMONIC PERTURBATIONS FOR MOM 
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CENTIMETERS 


COEFF. MODEL: 13x10“ 5 /^ 2 



Figure 2. Harmonic Perturbation by Degree for Radial, Along, and Cross Velocity 
Components for MO 
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MARS 50 x 50 HARMONIC MODEL 
MO: 70 DAYS DSN DOPPER 60 SEC. 0.03cm/sec. 



2 10 20 30 40 50 

DEGREE 

Figure 3. Error Spectra Per Degree 
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ON SEASONAL VARIATIONS OF MARS' GRAVITATIONAL FIELD 


B. Fong Chao 
David Parry Rubincam 


OBJECTIVE 

The objective of this work is to study, based on simplified models, 
the variation in Mars' gravitational field caused by the seasonal 
exchange of CO 2 mass between the Martian atmosphere and polar caps. 


BACKGROUND 

It is well-known that Mars has seasons like the Earth. However, the 
colder atmospheric temperature of Mars (some 80°K lower than the 
Earth's) and its predominantly CO 2 atmosphere (more than 95%) con- 
spire to create a large exchange of CO 2 mass between the atmosphere 
and polar caps. The mass exchange occurs in seasonal cycles: a 
great quantity of CO 2 condenses to form (or add to) the polar caps 
in winter and sublimes into the atmosphere in summer. The mass 
involved is at least 8 x 10 15 kg, which is about 25% of the total 
mass of the Martian atmosphere, or equivalent to a cap of solid CO 2 
extending to 45° latitude changing its thickness by some 20 cm. In 
contrast, the Earth's atmosphere only varies seasonally by about 1 
x 10 kg (or about 0 . 02 % of the total), primarily due to variation 
of water vapor content. 

The CO 2 mass redistribution inevitably changes Mars' gravitational 
field. As an order-of -magnitude assessment, for the gravitational 
effect of any mass redistribution on a planetary scale, the deciding 
factor is the ratio of the net redistributed mass to the total mass 
of the planet. For the Martian C 02 ~exchange effect, the ratio is 
8x1015/6 .4x10^3 t or about 10 " 8 . No geological or geophysical 
phenomenon on Earth is known to produce a relative mass redistribu- 
tion of nearly that size in a year. The post-glacial rebound of 
the mantle, the strongest southern oscillation/El Nino events in 
the atmosphere and the ocean, and the greatest earthquakes are all 
at least 2 or 3 orders of magnitude smaller. Since even these 
geophysical phenomena are already being observed gravitationally by 
Earth satellites (e.g., Rubincam, 1984; Gross and Chao, 1985), we 
can expect the C 02 -exchange effect to be detectable from future 
Mars orbiter missions. If so, the observed gravitational field 
changes can provide gross constraints on the meteorological models 
for Mars. 
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RECENT ACCOMPLISHMENTS 


Suppose the planet's surface is subject to a time-varying change in 
the surface density Ao(fl,t). The resultant change in the gravita- 
tional field, expressed in terms of its harmonic coefficients (or 
the "Stokes' coefficient") of degree i and order m can be shown to 

ho 


« c i.* ,s iJ (t > ■-^r^ (a ' t) vj dQ (1) 


where R is the radius, M the mass of the planet, Y im is the nor- 
malized spherical harmonic function, and the integration is oyer 
the entire surface (with surface element dn) in an Eulerian 
description. For convenience we shall also translate these changes 
into equivalent rms changes in the geoid height N using Bruns' 
formula 


4N, m (t) = R | 4(C ln +1S )m )(t) |. (2) 


AN^mU) in millimeters provides a comfortable measure of the effect 
of Ao(ft,t) on the geoid. It can also be shown that the elastic 
yielding of Mars is small so that the solid Mars can be treated as 
a rigid body. 

We follow the seasonal journey of the CO 2 mass and consider sepa- 
rately the contributions of two phases to equation (1): (i) the 
waxing and waning of the solid CO 2 in polar caps, and (ii) the 
geographical distribution of the gaseous CO 2 in the atmosphere. 
The two phenomena, of course, accompany each other inasmuch as the 
total CO 2 mass is constant. The total (^-exchange effect is just 
the sum of the two contributions. 


(i) We assume that the seasonal condensation and sublimation are 
uniform on the polar caps, and that the caps are circular and 
symmetrical with respect to the poles. Equation (1) then gives 


AC 10 (t) = 


1 


\/21+l (1-cos/c) 


1 AM (t) 

PjOO d/i 


COS/C 


M 


* (-1) 1 

same 

. n 

. 


(3) 


where P £ is the ordinary Legendre function of degree z , <=45° is 
the angular extent of the polar cap, AM is the total COjj mass 
accumulated on the cap at any given moment, and the subscripts n 
and s indicate the northern and southern hemisphere, respectively. 
Here we shall only consider a few of the lowest degree terms by 
evaluating their maximum amplitudes which, occur at the peak of 
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formation of the southern cap (t=t 0 ). Thus substituting into 
equation (3) the value for AM, AM(t Q ) = 8xl0 15 kg, and R=3390 km, 
the result shows that the largest cnange is in J 2 , which sees a 
maximum relative seasonal change of 3.8 ppm (parts per million). 
In comparison, the observed secular change in the Earth's J 2 (due 
to post-glacial rebound) is 0.024 ppm per year (Rubincam, 1984), 
which is two orders of magnitude smaller. The results are sum- 
marized in Table 1. 


(ii) The Martian surface topography will influence the way the 
changing atmospheric CO 2 mass distributes itself geographically; 
that, in turn, will change the gravitational field. Its effect on 
the gravitational field can be shown to be 


« c i. +,s i^w ■ - 


2-5. 


mO 


21+1 


RH 1m AM(t) 
Z M 


( 4 ) 


where H Am is the "effective topography" defined as the departure of 
the true topography from the geoid; Z is the scale height of the 
Martian atmosphere. Again, the largest effect is on J 2 . The 
results are summarized in Table 2. 


SIGNIFICANCE AND FUTURE EMPHASIS 

The variations in Mars' gravitational field due to the CO 9 mass 
exchange are generally very large compared with their terrestrial 
counterparts. Nevertheless, whether they can be observed by the 
upcoming MO (Mars Observer mission, due to be launched in 1990) is 
presently uncertain. A preliminary simulation (F. Lerch, private 
communication, 1985) has indicated that, using Doppler tracking data 
alone, the change in J 2 can be marginally detected, provided 
that all favorable conditions prevail throughout the MO lifetime 
(1 Martian year). However, if in addition the surface altimetry 
data (as differences at crossovers) are available, then the accuracy 
in the gravitational field determination is expected to be greatly 
improved. It may then be possible to detect changes not only in J 2 
but in other components as well. A more realistic and more detailed 
simulation is currently underway (see Lerch, this issue). If the 
above analysis is confirmed, the observed changes in the gravita- 
tional parameters can place constraints on the meteorological 
models of Mars. 
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Table 1 . Maximum Gravitational Effect of the Changing Polar Cap C0 2 (equation 3) 


t 

ACfo( X 10 -8 ) 

AJ f (Xl0 8 ) 

ANfo(mm) 

Remarks 

2 

0.34 

-0.75 

11 

Largest term. 

3 

-0.15 

0.40 

5.0 

Progressively smaller. 

4 

0.031 

-0.094 

1.1 

1 

5 

0.030 

-0.10 

1.0 



Table 2. Maximum Gravitational Effect of the Changing Atmospheric C0 2 (equation 4) 


t 

m 

AC fm (xl08) 

AS fm (xI0 8 ) 

AJ,(xl0 8 ) 

AN fm (rnrn) Remarks 

1 

0 

-0.053 

— 

0.092 

18 I 

i 

| Corresponding to 2.6 

> mm of shift in center 

1 

1 

-0.0023 

0.056 

— 

1.9 ) 

\ of mass 

2 

0 

-0.080 

— 

0.18 

2.7 

Largest term, about 
1/4 that of Table 1, 
but out of phase 

2 

1 

0.012 

-0.018 

— 

0.72 

Corresponding to 50 
cm of polar motion 
excitation 

3 

0 

0.0051 

— 

-0.014 

0.17 

Small compared with 
Table 1 

2 

2 

-0.034 

-0.028 



1.5 1 

1 Large because of 

3 

3 

0.0093 

-0.027 

— 

1.0 | 

( Tharsis 
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CONSTRAINTS ON VENUS' LITHOSPHERE RHEOLOGY FROM TECTONIC FEATURES 


Maria T. Zuber 


OBJECTIVE 

The geometries of tectonic surface features provide a direct 
indication of the rheological structure of the lithosphere at the 
time of deformation. Recent radar images of the surface of Venus 
reveal a variety of features which are interpreted to have formed 
in response to tectonic stresses. In this study, characteristics 
of these features are used to place constraints on theoretical 
models of the mechanical structure of the Venus lithosphere. 


BACKGROUND 

Radar images obtained from Pioneer Venus (Masursky et al . , 1980, 
Pettengill et al., 1980) earth-based observations (Campbell et al., 
1983, 1984), and the Venera 15/16 spacecrafts (Barsukov et al., 
1986; Basilevsky et al., 1983) indicate that the Venus surface has 
undergone extensive tectonic activity. In particular, the highlands 
regions of Venus (radius > 6053 km; Masursky et al., 1980) contain 
numerous examples of both extensional and compressional deformation 
(Head, 1986). A conspicuous extensional feature, illustrated in 
Figure 1, is the Beta Regio Rift zone. As observed in high resolu- 
tion Arecibo radar images (Campbell et al., 1984), the rift con- 
sists of a central depression and uplifted flanks with an average 
flank-to-f lank width of 100-200 km. Within the central depression 
are alternating radar-bright lineations which are interpreted as 
faults (Campbell et al., 1984) with widths which range from 10 to 
20 km. Figure 2 shows a group of prominent compressional features, 
termed ridge belts, which were revealed by Venera 15 and 16 imaging 
radar (Barsukov et al., 1986; Basilevsky et al . , 1986). The ridge 
belts are north-south trending linear ridges which have a regular 
spacing perpendicular to strike of approximately 300 km (Zuber, 
1986), within which are contained systems of sub-parallel trending 
ridges and grooves with widths of 10-20 km. The geometries of 
these features and others on the Venus surface illustrate that 
extensional and compressional deformation in many instances 
developed with two length scales, one of a few tens-of-km and 
another on the order of hundreds-of-km. 
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The parallel strikes and regular or characteristic spacings of 
tectonic features on Venus suggest that their length scales of 
deformation are controlled by a dominant wavelength arising from 
the growth of compressional and extensional instabilities in the 
lithosphere. Instability growth may occur in a horizontally- 
stressed medium which is stratified in viscosity (strength) or in 
which viscosity is continuous and decreases with depth (Biot, 1957; 
1960). Deformation arises under conditions for which initial 
perturbations along free surfaces and/or interfaces between regions 
of high and low viscosity amplify with time. The style of deforma- 
tion associated with instability growth is determined by the nature 
of the initial perturbations. In either an extending or compressing 
medium in which initial perturbations are distributed randomly, 
deformation develops periodically at the dominant wavelength, Xj 
(e.g., Fletcher, 1974; Smith, 1975). In an extending layered 
medium in which an initial perturbation is localized, deformation 
consists of a central depression and uplifted flanks, and is 
morphologically analogous to a rift zone (Zuber and Parmentier, 
1986). The dominant wavelength of the extensional instability 
controls the width of the rift. In both cases the dominant wave- 
length is a funcion of the mechanical structure of the medium, and 
in particular the thickness of the strong layer. Thus, instability 
growth may describe both the periodic spacings of extensional and 
compressional features and the width of an isolated rift. 

In a study of deformation due to extensional instability applied to 
the Basin and Range Province of the western U.S., Zuber et al . 
(1986) have shown that two scales of deformation may arise if the 
lithosphere consists of two high viscosity layers separated by a 
weaker layer. For the earth's continental lithosphere, the strong 
layers are believed to correspond to the upper crust and a region of 
strength in the upper mantle, while the weak layer is thought to 
correspond to the ductile lower crust. It is suggested that the 
existence of two scales of deformation on Venus indicates that its 
lithosphere has also differentiated into a crust that is chemically 
distinct and weaker than the underlying mantle. 

Figure 3 shows one possible strength stratification for the Venus 
lithosphere assuming that the thickness of the crust (shown by the 
dashed line) and the compositions of the crust and mantle are 
similar to those for terrestrial ocean basins. The regions in 
which strength increases linearly with depth correspond to the part 
of the lithosphere which deforms by brittle fracture. In other 
areas deformation occurs by ductile flow and strength decreases 
rapidly with increasing temperature (depth). The left and right 
sides of the diagram show the strength of the lithosphere, defined 
as the difference between the maximum and minimum principal effec- 
tive stresses, for compression and extension, respectively. Note 
that for both extension and compression the upper crust and mantle 
are stronger than the lower crust. If this scenario applies to 
Venus, then the spacings of small and large scale tectonic features 
may be controlled by the thicknesses of the strong regions of the 
upper crust and upper mantle, respectively. 
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RECENT ACCOMPLISHMENTS 


To test the hypothesis that tectonic features on Venus formed in 
response to multiple scales of extensional and compressional insta- 
bility, simple models of the Venus lithosphere consisting of two 
strong or high viscosity layers separated by a weaker layer have 
been developed. Results show that the lithosphere can be unstable 
in both extension and compression and that a range of conditions 
exist for which two scales of deformation will occur. Figure 4 
illustrates how the dominant wavelengths are controlled by the 
thickness of the strong crustal layer (hj) and some of the physical 
properties of the medium for both compression and extension. The 
two lines in each diagram correspond to the short and long wave- 
lengths of deformation. If the width or spacing of a tectonic 
feature is known, then these diagrams can be used to constrain the 
thickness of the strong crustal layer. Similar relationships can 
be established for other physical properties in the problem. 
Results show that the allowable range of strong crustal layer 
thicknesses is 3<hj<20 km and that total crustal thicknesses fall 
generally in the range 5<h c <30 km. A strong region of the mantle 
only a few km thick is required to produce a long wavelength 
instabi lity. 

In conclusion, compessional and extensional instability growth is a 
viable mechanism for the formation of many tectonic features on 
Venus. The results suggest that at least in some areas, the Venus 
lithosphere consists of two regions of strength, analogous to the 
strong upper crust and upper mantle of terrestrial continents. 


FUTURE EMPHASIS 

Previous analyses (e.g.. Smith, 1975) have shown that a compressing 
layered medium is always unstable (and therefore always deforms at 
the dominant wavelength) in compression, while in extension deforma- 
tion may be either unstable or stable (corresponding to uniform 
thinning of the medium). On Venus it may be possible to define a 
plausible range of conditions for which the lithosphere is stable 
in extension. This would place important constraints on the 
strength of the lithosphere and on the origin of the tectonic 
features. 
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Figure 1 . Topographic profiles and sketch map of bright and dark linear features mapped from 
Arecibo radar image of the Beta Regio Rift (from Campbell et al., 1984). The radar- 
bright lineations within the central depression of the rift are interpreted as faults. The 
width of the rift and the spacing of the lineations define two length scales of 
extensional deformation. 



Figure 2. Sketch map showing the distribution of ridge belts in part of the northern hemisphere 
of Venus (after Barsukov et al., 1986). Within the ridge belts, smaller scale ridges 
and grooves strike approximately parallel to the belts. The spacings of the ridge belts 
and the ridges and grooves define two length scales of compressional deformation. 
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Figure 3. Strength of the Venus lithosphere with depth for compression and extension assuming 
a crustal thickness and crust and mantle compositions similar to those of terrestrial 
ocean basin. The strength envelope was constructed for the geothermal gradient 
shown on the right. 
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normalized to the thickness of the crustal layer hj. The two lines correspond to short and long wavelength instabilities. 




CHAPTER V 


GEOPHYSICAL TECHNIQUES DEVELOPMENT 


OVERVIEW 

This chapter discusses a variety of subjects concerned with 
measurement system development, computer software development and 
enhancement, and the development and evaluation of new measurement 
techniques. Specifically, a detailed discussion is made of the 
modification, through use of new technology, of the Very Long 
Baseline Interferometer (VLBI) system, supporting computer software 
development, and mathematical modeling based on the analysis of VLBI 
data. A discussion is provided of the calibration procedure used 
on new Satellite Laser Ranging (SLR) Systems prior to their intro- 
duction into the global SLR network used to track the Lageos 
satellite. Specialized techniques such as computer imaging techni- 
ques as applied to ocean topography and geophysical data are 
described. Analyses to evaluate the capabilities of new measure- 
ment techniques such as (1) the direct measurement of the high 
frequency components of the earth's gravity field through use of a 
spaceborne gravity gradiometer and (2) the global measurement of 
crustal deformation, land topography, etc., through the use of a 
spaceborne geodynamics laser ranging system are both discussed in 
this chapter. Finally, a discussion of research on soil tempera- 
ture is included. 

Contributors to this chapter are: Thomas A. Clark, Steven C. Cohen, 

Theodore L. Fel sentreger, John M. Gipson, M.W. Hayes, Edwin Himwich, 
Werner D. Kahn, Ronald Kolenkiewicz, Goran L. Lundqvist, Chopo Ma, 
Carole Nethery, Barbara H. Putney, James Ray, and Jean E. Welker. 
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MICROWAVE HOLOGRAPHY 


John M. Gipson 


OBJECTIVE 

During 1984, Very Long Baseline Interferometry (VLBI) was used to 
measure the complex (amplitude and phase) far field pattern of 
three antennas used by the Crustal Dynamics Project with the aim of 
using microwave holography to obtain information about these 
antennas . 


BACKGROUND 

Microwave holography is a technique for obtaining information about 
an antenna from its far field pattern. One can find the far field 
pattern by modeling the surface of the antenna as a collection of 
elementary current dipoles. A dipole is a simple kind of antenna 
whose radiation pattern is known. The strength of the radition 
pattern is proportional to the dipole current. By adding up the 
radiation pattern from all of the elementary dipoles one can obtain 
the far field of the antenna. In microwave holography one inverts 
this relation to obtain the effective dipole distribution from the 
far field pattern. This information can be presented in the form 
of maps which display the strength of the signal on the surface of 
the antenna (Fig. 1), and deformations of the antenna surface (Fig. 
2). One can use this information to improve the performance of an 
antenna. Ideally, the surface of an antenna should be a perfect 
parabola. Deviations of the surface from the ideal result in a 
loss of power while transmitting, and a loss of sensitivity while 
receiving. An rms deviation of 1/16 of the transmitting (observing) 
wavelength results in a 50% loss in power (sensitivity). By 
measuring the surface of an antenna one can determine which areas 
need adjustment, and by how much. 

The measurement of the far field is done with respect to a reference 
antenna. All-in-all, four experiments were done. Table 1 sum- 
marizes these experiments. 

During the experiments the reference antenna tracked a celestial 
source. The antenna to be measured pointed at the celestial source, 
and at an N-by-N grid of points centered on the source. Both 
antennas measure and record the strength of the received signal 
across a broad band of frequencies. These data are subsequently 
manipulated to obtain the far field of the antenna. This informa- 
tion in turn is used to obtain maps of the antenna. 
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Table 1. 



Measured 

Reference 

Date of Experiment 

Antenna 

Antenna 

January 12, 1984 

Mojave 

Goldstone Venus 

August 31, 1984 

Gilmore Creek 
(Fai rbanks) 

Mojave 

October 21, 1984 

Mojave 

Owens Valley 

October 21, 1984 

Owens Valley 

Mojave 

The resolution length R 

(the size of the 

minimum observable 


features) is determined by the antenna size 0, the sampling factor 
f, and by the grid size N: 

R = D/(N * f). 


Here f is a number between 0 and 1 which determines how much of 
the final map the antenna occupies. If f = 1, the antenna just 
fills the map. In order to clearly define the antenna, a value of 
f = 0.5 was used. (See figures 1 and 2). In retrospect, this 
value is excessively conservative. It appears that f = 0.8 would 
be a good tradeoff between observing efficiency (number of points 
required to achieve a given image resolution) while avoiding edge 
effects. An 11 by 11 grid was used in the first experiment, 
and a 13 by 13 grid in the remaining experiments. The angular 
spacing a between the grid points is relatd to the antenna size D, 
the observing wavelength A, and the sampling factor: 


A = f A/D. 


RECENT ACCOMPLISHMENTS 

The analysis of the four experiments summarized in Table 1 is now 
complete and will come out shortly as a Goddard Technical 
Memorandum. 

Preliminary analysis of the data was done at Haystack, and the 
results transmitted to Goddard. The far field pattern of the 

antenna was extracted from the data at Goddard using FRNGE, a 
program written by Alan Robers of Haystack. The far field pattern 
was then manipulated to obtain the effective dipole current 
distribution. 
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The amplitude of the current distribution is proportional to the 
strength of the signal on the antenna. By measuring it, one can 
directly examine the effects of feed and support leg shadowing, 
areas of low reflectivity, and the aperture illumination function. 
Fig. 1 shows a plot of the power distribution for the antenna at 
Gilmore Creek, in Fairbanks, Alaska. Generally, the power peaks in 
the middle of the antenna and falls off smoothly. Feed box shadowing 
causes a 4 db drop in the signal strength at the center of the 

antenna. The legs which support the feed run in the NS and EW 

directions, and the effects of support leg shadowing are clearly 

visible. Prior to this experiment, a panel in the SW quadrant of 

the antenna was removed. Its absence caused a 2.6 dB loss in this 
area. 

The phase of the current distribution yields information about the 
reflector surface profile and the feed location. In the case of a 
reflector antenna, the phase of the current distribution is 
proportional to the distance from the antenna surface to the feed. 
Deviations of the phase from its expected value indicate areas of 
the antenna where the distance from surface to feed differs from 
the ideal case. Such deviations can be due to two effects: a 

deformed surface, or a misaligned feed. 

Surface deformations due to badly formed or aligned panels show up 
as local deviations in the phase. Large scale surface deformations 
due to thermal effects or gravitational loading lead to large scale 
deviations in the phase. Figure 2 shows a map of the surface profile 
of the Gilmore Creek antenna. The antenna is smooth with a few 
exceptions. There are sharp gradients at 1 o'clock and at 4 o'clock. 
There is another sharp gradient at the 7 o'clock position. This is 
probably connected with the missing panel in this area. The power 
of the antenna could be improved by about 8% by removing the slight 
surface deformations. 

Feed alignment errors lead to large scale deviations with 
characteristic signatures. Lateral feed misalignment introduces a 
phase gradient across the antenna: The feed is closer to one side 

of the antenna than the other, and hence signals reach one side 
before they reach the other. A radial feed alignment error causes 
a phase gradient from the inside of the antenna out. In this case, 
the feed is closer to (or further from) the inner regions of the 
antenna than it should be with respect to the outer regions. The 
characteristic effects of feed alignment errors makes it possible 
to detect, and compensate for them in the analysis of the phase 
maps. None of the antennas measured showed significant feed 
alignment errors. 
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SIGNIFICANCE 


These experiments demonstrated the feasibility of using VLBI to 
perform microwave holography. In principle, one could use the 
results of these experiments to adjust the panels and the feeds of 
the measured antennas to improve their performance. The expected 
improvements in efficiency range from 5% to 28% at X-band. (The 
VLBI sensitivity goes as the square root of the single antenna 
sensitivity, so a great deal of sensitivity is not lost). At higher 
frequencies the improvement would be more dramatic. 


FUTURE EMPHASIS 

The analysis of these experiments suggests two other experiments. 
It would be interesting to make higher resolution maps. It would 
also be interesting to study the effects of gravity by making maps 
with the antenna pointing at a number of different directions. 

The effects of individual panels are not clearly apparent in any of 
the maps of this report. For example, in-the Gilmore Creek antenna, 
the resolution length is twice the length of a panel. If f were 
increased to 0.8 and the grid size to 16 the resolution length 
would be comparable to the panel size. The increased grid size 
presents problems with the signal to noise ratio (SNR), and with 
the time required to acquire the data. For the experiments in this 
report, the SNR of the on-source points was typically 500, while 
the SNR of the points at the edge of the grid was about 20. As the 
grid size increases, the points at the edge are increasingly far 
away from the beam center. Therefore, to achieve a moderate SNR of 
10 one would have to sample for a longer period of time. This 
increases the time required to acquire data. It may no longer be 
possible to complete an experiment in a single apparition of a 
celestial source. However, the data from two or more nights could 
be combined to make a map with the desired resolution, but with 
possible contamination from thermal effects. 

The second suggested experiment is in some sense the converse of 
the first. It would be interesting to take 'snapshots' of the 
antenna at several different positions to study the effects of 
gravity. By sampling at f = 0.8 on an 8 by 8 grid we obtain a map 
of the same resolution as those reported here. This can be done in 
under an hour. One could therefore take several maps during a 
single apparition of a celestial source. 
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Exp 441. Gilmore Creek. Power Plot. 

Figure 1 . The center of the antenna is marked with an x. The physical region of the antenna is 
indicated by the circle. The power is in dB. 
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Exp 441. Gilmore Creek. Surface Plot. 


Figure 2. The center of the antenna is marked with an x. The physical region of the antenna is 
indicated by the circle. The surface deviation is in millimeters. 
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USE OF PHASE DELAY OBSERVATIONS FOR GEODETIC DETERMINATIONS 


Jim Ray 


OBJECTIVE 

Normally geodetic VLBI measurements are made using the group delay 
observable, the fundamental precision of which is limited by the 
radio bandwidth that can be synthesized (i.e., spanned using a set 
of relatively narrow frequency channels). This report summarizes 
work at GSFC to greatly improve VLBI measurement precision by using 
the phase delay observable in addition to the group delay. 


BACKGROUND 

The VLBI geodetic technique uses radio emissions from celestial 
sources observed simultaneously from two remote sites to determine 
the phase difference between the relative wavefront arrivals at the 
two locations. In practice, the phase determination is made rela- 
tive to stable frequency standards maintained at the sites. The 
phase difference is a direct and very precise measure of the base- 
line separating the observing stations projected into the direction 
to the source. As this measurement can only be made modulo an 
integral number of wavelengths, interpretation of the observation 
is inherently ambiguous by an unknown number of wavelengths. If 
however, phase measurements are made simultaneously at several 
different frequencies, then the derivative of the phase with respect 
to frequency will be free of cycle ambiguity problems. This latter 
data type is known as the group delay. Currently, VLBI geodesy 
relies almost exclusively upon group delay measurements. 

While the frequency derivative of the phase does allow the cycle 
ambiguity problem to be circumvented, the measurement precision 
degrades by the ratio of the bandwidth spanned in forming the group 
delay to the frequency of the band center. For the frequency 
sequences used by the Crustal Dynamics Project (CDP), the group 
delay observations are less precise than phase delays by a factor 
of more than 20. 
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Observations have been scheduled so as to dwell within each cluster 
for about 2 hours at a time. This approach of minimizing large 
angular sweeps should enhance the chance for maintaining phase 
connection throughout the experiment as well as allow a study of 
relative phase stability as a function of source-source separation 
angle. 

The CDP experiment of March 8, 1985 produced resolvable S- and X- 
band phase delays. This was accomplished in a two-step process 
that first relied upon a very good group delay solution to determine 
values for the relative site positions, clock behavior, zenith 
tropospheric delay, and positions for some of the sources. Using 
the adjusted values for these parameters, resolution of the phase 
ambiguities was relatively straightforward. The phase delays then 
gave a solution having a postfit rms residual of 14 ps (4 mm). 
Comparison of the phase solution with the group solution from the 
same experiment shows agreement in baseline length within 2 mm. 
The formal uncertainty on the baseline length from the phase solu- 
tion is 1 mm (1 sigma). There is, however, a significant dif- 
ference in the orientation of the two baseline determinations. The 
reason for this appears to be a low-level phase calibration error 
in the OVRO receiver system. The problem had not been noticed in 
the usual group delay solutions, but became evident through close 
inspection of the phase delay residuals. Both the phase and group 
solutions from this experiment are consistent with the average 
result of a group of five mobile experiments also conducted during 
March 1985, but the formal uncertainty for the mobile experiments 
is substantially larger. 


SIGNIFICANCE 

Thus far, the high accuracy of the phase delay measurements has 
helped to pinpoint non-random variations in instrumental calibration 
systems that were not detectable using group delays. For the same 
reason, however, we have not yet established the level of atmos- 
pheric phase noise for intermediate length baselines. Apparently, 
it is low enough, at least some of the time, to permit phase connec- 
tion over distances of at least 250 km. With continued improvement 
in instrumental stability, it appears that measurement of 
California baseline lengths to about the millimeter level should be 
feasible using phase delays. 


FUTURE EMPHASIS 

A current goal of this effort is to extend this technique to a 
three-baseline network by including the Hatcreek Observatory 
together with Mojave and OVRO. The longest baseline in this 
example is 729 km. In addition, an attempt will be made to phase 
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With recent improvements in the CDP receiver technology, frequency 
standards, calibration techniques, and a priori site and source 
positions, group delay solutions for many baselines of intermediate 
length have achieved postfit residual values (rms) of one-quarter X- 
band wavelength or less (i.e., about 30 ps or 1 cm). In these 
cases, it is reasonable to expect a group delay solution adjustment 
for variable parameter values (e.g., site positions, zenith tropos- 
phere delay, site clock differences) should be adequate to permit 
an unambiguous resolution of the phase delay ambiguities. In this 
bootstrapping approach, the full potential of the precise but 
ambiguous phase delays may be realized. 

While interferometric phase observables are already used routinely 
for astronomical observations, those measurements are normally made 
possible by distributing timing signals to all stations from a 
central facility. The physical links required for such phase- 
coherent networks necessarily limit their scale to a few tens-of- 
kilometers, much too small for most applications of interest in the 
CDP. Demonstration that phase delay observations could be success- 
fully utilized in a VLB I configuration with independent clocks was 
made in 1978. Results for the 1.24 km baseline between the 
Haystack and Westford telescopes in Massachusetts showed sub- 
centimeter repeatability (Rogers et al., 1978). 


RECENT ACCOMPLISHMENTS AND RESULTS 

Recent efforts at 6SFC have attempted to extend the baseline scale 
over which phase delays can be utilized and to understand the 
phenomena which limit the application of the technique to still 
longer baselines. The first step in this effort was a series of 
experiments to measure the 12.6 km baseline between the Mojave Base 
Station and the DSS-13 antenna located in the DSN Goldstone complex. 
Consistency at the sub-centimeter level was achieved in the two 
horizontal components of the baseline while a roughly 2 cm result 
was obtained for the vertical component. The poorer quality of the 
latter component can be attributed to error in the tropospheric 
calibration, an effect largely independent of the data type used. 

The next step has been to attempt to phase connect the Mojave 
antenna to the 130 foot dish at the Owens Valley Radio Observatory 
(0VR0), a separation of about 245 km. This much longer baseline 
represents a formidable challenge to the technique, mostly on 
account of the greater phase fluctuations due to more highly 
uncorrelated atmospheric columns over sites widely separated from 
one another. In an attempt to minimize problems of this type, an 
innovative schedule has been developed using tight clusters of 
sources each having internal separations of less than 20 degrees. 


163 


connect Hatcreek to a mobile station based at Quincy. While this 
baseline is considerably shorter than the Mojave-OVRO separation 
(100 km verus 245 km), the attempt is complicated by the relative 
insensitivity of the mobile receiving systems compared with the 
fixed stations. A highly precise determination of the Hatcreek- 
Quincy baseline is desirable in order to map baseline measurements 
from one site to the other for the purpose of better estimating 
tectonic motions between northern and southern California. 


REFERENCE 

Rogers, A.E.E., C.A. Knight, H.F. Hinteregger, A.R. Whitney, C.C. 
Counselman III, I. I. Shapiro, S.A. Gourevitch, and T.A. Clark, 
"Geodesy by Radio Interferometry: Determination of a 1.24-km Base 

Line Vector with 5-mm Repeatability," J. Geophys. Res ., 83 , 325- 
334, 1978. 
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GLOBL : GLOBAL LEAST-SQUARES SOLUTION SOFTWARE FOR VLBI DATA 


Edwin Himwich 


OBJECTIVE 

Since 1978, the VLBI Data Analysis Team of the Geodynamics Branch 
has had a continuing task to develop a software system that could 
be used to analyze the ensemble of geodetic VLBI data to recover 
estimates of geophysically and astrometri cal ly interesting quanti- 
ties in an efficient and flexible manner. From the standpoint of 
the Crustal Dynamics Poject the most important quantities to be 
extracted are the baseline lengths, that is, the distances between 
pairs of VLBI stations. 


BACKGROUND 

The program used for the operational analysis of single VLBI 

experiments, SOLVE, was written in 1976. A single experiment 
consists of about 24 hours of data. In order to get good estimates 
of station positions, SOLVE allows data analysts to study the 
behavior of the clocks and atmosphere during the experiment. 

To analyze data from two or more experiments, a new program, S0LV2 
was written in 1981. S0LV2 proved to be very reliable and is 

extremely useful for comparison with SOLVE for the purpose of 
verification. The program provided tremendous flexibility in its 
operation, but at the cost of being extremely cumbersome. It 

eventually became clear that in order to meet the analysis needs of 
the project, a new program would have to be developed. 


RECENT ACCOMPLISHMENTS AND RESULTS 

The program GLOBL has been written to make least squares parameter 
estimates from several VLBI experiments. It uses the best features 
of S0LV2 and other programs based on SOLVE to provide a much more 
efficient and much less cumbersome package than S0LV2. It runs 
much faster than S0LV2 and requires much less setup. 

Figures 1 and 2 show the results from a single GLOBL run to estimate 
nutation errors. Figure 1 shows daily offset of nutation in 
longitude; figure 2 shows offset in obliquity. The data used for 
this solution includes 167 experiments, all of the IRIS and CDP 
data for 1984 and most of 1985. The most interesting feature of 
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this data is the apparent annual trend. There are hints of other 
effects in the data. In particular, there appears to be some small 
scale structure with a period of about 14 days. A slope appears to 
be present in the longitude data as well. More analysis will be 
necessary before precise estimates of these effects can be made. 

Some of the more important features of GLOBL are: 

A. Saving the Data in a High-Speed Access form (Superfiles) 

VLBI data is normally stored in what is known as database 
form. The database system provides a high level of accountability 
for the steps taken in processing the data. This self-documenting 
storage form has a relative low access rate. In GLOBL, advantage 
is taken of the fact that data is not included in a global solution 
until its processing has reached a stable point. GLOBL never 
modifies the data, so the need to update the permanent documentation 
is avoided. This allows the data to be stored in a higher speed 
access form, known as Superfile form. 

B. Arc-Parameter Elimination 

GLOBL, like S0LV2 before it, takes advantage of the large 
scale structure of the normal equations! . A global solution of 
about half the Mark III data set would use data from about 200 
single day experiments. The direct solution of the normal equations 
for this case would require the inversion of an approximately 10,000 
by 10,000 element matrix. Arc-parameter elimination notes that 
most of this matrix is zero and reduces the matrix inversions to 
one 120-by-12U and 400 50-by-5O matrices. The time required for 
matrix inversion goes up roughly as the cube of n for an n-by-n 
matrix, all other things, such as everything being able to fit in 
core, being equal. Thus arc-parameter elimination reduces the time 
spent in matrix inversion processes by a factor of almost 20,000. 
Without arc-parameter elimination, 2 years would be required for 
just matrix inversions. 

C. Statistics 

In order to understand the effect of solving for different 
parameters or of treating parameters differently, it is useful to 
have statistics available from the solution. The scaled formal 
error for each parameter and also residuals by baseline and source 
are provided. The residual statistics are calculated individually 
for each experiment as well as for the overall solution. 

D. Command File 

GLOBL solutions are controlled by a command file. The command 
file contains all the information that is required to run a solution. 
This is a fairly small amount of information. In addition to about 
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two dozen control parameters, the name of each experiment to be pro- 
cessed must be provided. This compares to S0LV2 which required 
almost two dozen control parameters to be specified for each experi- 
ment and required a program to prompt for the values of the 
parameters. The simplification in GLOBL which reduces the setup 
burden is to break parameters that can be solved for into different 
classes. GLOBL than provides control over how each class will be 
treated. There is a tremendous loss in terms of absolute control 
of the solution. However, in practice most of that control is not 
needed. The resulting savings in setup time allows more time to be 
spent on more important activities, such as interpretation. 

E. Archiving of Results 

GLOBL makes use of SOLVE's solution archiving system. This 
allows the results of solutions to be saved for use with the 
automatic report generator. The report generator will make tables 
documenting the time evolution of baseline lengths and earth 
orientation parameters. This will be expanded to include tables of 
other parameters that can be estimated by GLOBL. 


FUTURE EMPHASIS 

The goals for the immediate future are to expand the types of 
parameters that can be solved for and to improve the running speed 
of the program. The parameters to be solved for can be easily 
increased to include the coefficients of the nutation and solid- 
earth tide models. This will allow direct solution for adjustments 
to both models. A decrease in running time of approximately 33% is 
anticipated with a straight-forward streamlining of the matrix 
operati ons. 
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Figure 1 . Corrections to Nutation in Longitude from Mark in VLBI Data vs Date 
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Figure 2. Corrections to Nutation in Obliquity from Mark m VLBI Data vs Date 



CREATION OF A GLOBAL GEODETIC NETWORK USING MARK III VLBI 


Chopo Ma 


OBJECTIVE 

For many geodynamic studies it is necessary to have a well-defined, 
well -controlled terrestrial reference frame of known precision. 
VLBI is a means for providing such a reference frame. 


BACKGROUND 

The Mark III VLBI system was developed by NASA to support research 
in geodynamics and has been used in high precision geodesy since 
1979 (Clark, 1979; Rogers et al., 1984; Carter et al . , 1985; Clark 
et al., 1985). By the end of 1984 Mark III data had been acquired 
from 15 stations which are either dedicated to or available for 
geodetic measurements. These are listed in table 1, The Haystack 
Observatory is the site of two stations, the 37-m Haystack antenna 
and the 18-m Westford antenna. Only Harvard, Haystack, Onsala, 
Owens Valley and Westford were used before 1983. In addition, five 
other stations have contributed data on an occasional or one-time 
basis: Chilbolton, England; Effelsberg, FRG; Maryland Point, 

Maryland; National Radio Astronomy Observatory, Western Virginia; 
and Robledo, Spain. Five stations regularly participate in IRIS 
measurements (Harvard, Richmond, Westford and Wettzell every 5 
days; Onsala monthly). Except for Richmond, these stations also 
participate with the others in observing sessions of the Crustal 
Dynamics Project. 

The status of the stations varies. Kashima, Onsala and Wettzell 
have permanent VLBI staff and Mark III equipment and are supported 
by national research agencies. Harvard, Haystack, Mojave, Richmond 
and Westford likewise have full-time staff and Mark III systems and 
are funded largely or entirely by NASA and NGS. At Algonquin, 
Gilmore Creek, Hat Creek, Kauai, Kwajalein, Owens Valley and 
Vandenberg, Mark III equipment and personnel are made available 
when required by the Crustal Dynamics Project observing program. 
Gilmore Creek and Vandenberg are operated solely for the Project 
while the remainder primarily support other unrelated activites. 

The 15 stations form a network whose relative positions are entirely 
determined by Mark III VLBI data. Since VLBI has no direct tie to 
the center-of-mass or to the spin axis, the origin and orientation 
of the network must be defined. The origin is defined by the 
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adopted Cartesian coordinates of the intersection of axes of the 
Haystack 37-m antenna. The orientation is defined by the subnetwork 
in use on October 17, 1980 (Chi Ibolton, Harvard, Haystack, Onsala 
and Owens Valley) and the a priori earth orientation parameters 
used to calculate the theoretical VLB I delays and rates for this 
observing session. These a priori values were quadratical ly 
interpolated from the BIH Circular D to each observation epoch, and 
the UT1 values were additionally modified to restore the short- 
period tidal variations. Any errors in Circular D for this period 
will be directly propagated into the solution and into the overall 
orientation of the network. The relative orientations of the 

baselines are not affected by errors in Circular 0 and would not be 
altered by a change in the reference day. Likewise, the baselines 
are not affected by a change in origin. 

To add new stations properly to the original network, it is 
necessary to observe simultaneously with at least four stations of 
which at least three must be in the original network. From the 
three or more old stations the complete earth orientation parameters 
can be estimated simultaneously with the coordinates of the new 
station(s) in the system of the original network. Thus, for 
example, Westford was added using the network of Harvard, Haystack, 
Onsala, Owens Valley and Westford from June 20, 1982. It is desir- 
able to link as many stations as possible through direct simul- 
taneous observations, but practical considerations limit current 
networks to seven stations. In addition, the requirement of mutual 
visibility of sources constrains the size of networks. Some of the 
networks used in 1984 were designated Polar (Haystack, Mojave, 
Gilmore Creek, Kashima, Wettzell), East Pacific (Mojave, Vandenberg, 
Gilmore Creek, Kauai, Kwajalein), West Pacific (Mojave, Gilmore 
Creek, Kauai, Kwajalein, Kashima), and North American Plate 
Stability (Westford, Algonquin, Harvard, Gilmore Creek). The first 
three networks were observed twice in 1984, each session lasting 
30-56 hours, and all will be repeated regularly for the duration of 
the Crustal Dynamics Project. 


RECENT ACCOMPLISHMENTS AND RESULTS 

Data from 16 stations acquired during 1980-84 by the Crustal 
Dynamics Project and POLARIS/IRIS were combined in a single, sequen- 
tially incremented least squares adjustment including 46070 observa- 
tions, each observation comprising a delay and a delay rate. Most 
individual data sets spanned 1 day. All data had ionospheric and 
tropospheric calibration. There were 132 global parameters 
estimated: all station positions except Haystack and all 44 source 
positions except the right ascension of 3C273B, whose adopted value 
defines the right ascension origin. There were 2499 arc parameters 
including earth orientation parameters from each session except the 
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reference day, polynomials and sinusoids to model station clock 
behavior, and zenith tropospheric delays. The models used 
generally followed the MERIT standards except for the absence of 
the K1 solid earth tide correction and ocean loading. Details of 
the data set and the solution are given in Ryan and Ma (1985). The 
overall fit of the solution was 120 ps in delay and 0.09 ps/s in 
delay rate. 

The distribution of observations by station is given in table 2. 
The total is twice the number of observations since each observation 
counts toward two stations. The observations from the POLARIS 
stations are most numerous, but even the limited use of Algonquin 
ties it to the network with 416 observations. The theoreticals do 
not include any a priori station drift because of tectonic plate 
motion, but the baselines for which data are available over the 
entire time span are not changing more than 2 cm/yr. The data from 
the Pacific stations come only from 1984. 


SIGNIFICANCE 

The adjusted source and station positions can be found in Ryan and 
Ma (1985). The positions of the 15 permanent VLBI stations are 
determined with 1-sigma formal statistical errors of less than 5 cm 
except for three stations in the Pacific. Based on internal and 
external tests of consistence, the true uncertainties are probably 
larger by a factor of 2. While the origins of the celestial and 
terrestrial coordinate systems are specified separately, the source 
and station positions are not independent. In particular, there is 
a correlation between the scale of baseline lengths and the 
declination of the sources. Except for some of the 1984 data from 
the Pacific, the polar component of baselines is generally much 
less than the component projected on the equatorial plane. The 
determination of source declinations, especially for sources near 
the equator, is correspondingly weak. 


FUTURE EMPHASIS 

The Crustal Dynamics Project acquires VLBI data from the global 
network on a regular basis. The quality of the data from more 
recent years indicates that the network could be tied together at 
the centimeter level. There are, however, improvements in the 
application of such models as global plate motions and nutation and 
in the calibration of the troposphere which must be made before 
this level of accuracy can be achieved. 
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Table 1. Global VLBI Network Stations 

Algonquin Radio Observatory, Ontario 

Gilmore Creek, Alaska 

Hat Creek Radio Observatory, California 

Harvard Radio Astronomy Station, Texas 

Haystack Observatory, Massachusetts 

Kashima Space Research Center, Japan 

Kokee Park STDN (Kauai), Hawaii 

Kwajalein, Marshall Islands 

Mojave, California 

Onsala Space Observatory, Sweden 

Owens Valley Radio Observatory, California 

Richmond USNO Timing Station, Florida 

Vandenberg, California 

Satelitenbeobachtungstation Wettzell, FRG 


Table 2. Distribution of Observations by Station 


Station 

Number of Observations 

Algonquin 

416 

Chilbolton 

1973 

Gilmore Creek 

3276 

Hat Creek 

708 

Harvard 

27745 

Haystack 

5975 

Kashima 

1819 

Kauai 

2110 

Kwajalein 

1185 

Mojave 

5686 

Onsala 

5956 

Owens Valley 

4507 

Richmond 

1743 

Vandenberg 

795 

Westford 

25560 

Wettzell 

2686 
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COLLOCATION RESULTS FOR THE NASA AND FRG LASER SYSTEMS 


Ronald Kolenkiewicz 


OBJECTIVE 

An objective of the NASA Crustal Dynamics Project is to obtain 
Satellite Laser Ranging (SLR) data from the Laser Geodynamics 
Satellite (LAGEOS). In order to achieve this goal new laser systems 
must be developed and compared to each other. This process of 
comparison is called collocation. The standard of measure for the 
NASA lasers is currently MOBLAS-7, located at the Goddard Space 
Flight Center in Greenbelt, Maryland. During the summer of 1985 
the newly developed laser system, MTLRS-1, owned by the Federal 
Republic of Germany (FRG), arrived in the USA to perform a colloca- 
tion with MOBLAS-7. The purpose was to assess the performance and 
precision of this new laser and to transfer the laser measuring 
standard to Europe. 


BACKGROUND 

In April 1985 the German satellite laser ranging system, MTLRS-1, 
was shipped to the NASA Goddard Space Flight Center in Greenbelt, 
MD. Here this system was tested by collocating the NASA mobile 
laser ranging system, MOBLAS-7. The object of a satellite laser 
ranging system collocation is to place two systems side-by-side, at 
a known orientation to each other, for the purpose of tracking the 
same target satellite and to compare the range data obtained. The 
two systems were precisely surveyed at a distance of 42.37 meters 
from each other. For the purpose of this experiment the MOBLAS-7 
laser was to be considered the standard by which the MTLRS-1 could 
be evaluated. Satellites to be tracked in order of their priority 
were LAGEOS, Starlette, and BE-3. The analysis presented here is 
for the LAGEOS satellite only. This satellite was selected since 
it is at an altitude which permits laser tracking data to be obtained 
for a continuous time period as long as 50 minutes, whereas the 
other lower satellites can only be tracked for a maximum of 15 
minutes. In addition, LAGEOS is going to be the prime target of 
MTLRS-1 when this system participates in the Mediterranean Campaign, 
conducted by the Working Group of European Geoscientists for the 
Establishment of Networks for European Research (WEGENER), beginning 
in 1986. 
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RECENT ACCOMPLISHMENTS 


During the period from May through July 1985, LAGEOS passes were 
tracked simultaneously at GORE by MOBLAS-7 and MTLRS-1. Using data 
from these passes as a starting point, a collocation analysis was 
conducted for the two systems. 

The collocation analysis technique uses the precision orbit 
determination computer program, GEODYN, to obtain a one pass (50 
minute maximum) LAGEOS orbit by using the full rate laser ranging 
data obtained solely by MOBLAS-7. The full rate data obtained 
during the same pass by MTLRS-1 is superimposed on the orbit, but 
in no way does it influence the MOBLAS-7 determined orbit. Range 
residuals to the orbit are thus obtained. Beginning at the first 
even minute before ranging data for the pass is obtained, 2 minute 
time intervals are chosen for the entire length of the pass. The 
mean of the range residuals and their standard deviations are 
calculated for MOBLAS-7 and MTLRS-1 during each of these 2 minute 
periods, and their value is considered to be at the time of the 
midpoint of the interval (the odd minute). In this manner 
corresponding 2-minute mean residuals from each system can be 
subtracted, where they exist, and the systematic trends between 
them may be evaluated. 

The first LAGEOS collocation pass candidate was taken on May 10, 
1985 at 0600 GMT by MOBLAS-7 and MTLRS-1 and the collocation 
concluded on July 16, 1985. During this time period 32 LAGEOS 
passes were tracked simultaneously by both systems. Of these, 13 
were selected as having sufficient data to be adequate for a 
collocation analysis. The selection criteria used to identify 
these passes follows: 

1. Each pass was to have at least five consecutive 2-minute 
mean points, simultaneously obtained from MOBLAS-7 and MTLRS-1. 

2. A valid 2-minute mean point could only be formed if there 
were at least 300 points from MOBLAS-7 and at least 30 points from 
MTLRS-1 during the 2-minute period. 

This selection criteria resulted in the further consideration of 
the 13 passes given in Table 1. None of the passes listed in this 
table were taken during the daylight hours, a condition that is 
always listed as a requirement for a successful collocation test. 
All of the passes are on the descending portion of the retrograde 
LAGEOS orbit, thus passing from northeast to southwest. It is 
always desirable to track ascending satellite orbits as well during 
a collocation of systems in order to aid in the determination of 
any azimuth or orientation related systematic errors that may exist. 
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Figure 1 shows LAGEOS collocation data for a typical pass taken on 
June 26, 1985 at 0400 hours GMT. The range residuals are plotted 
as a function of elapsed time from the beginning of the pass for 
both the M0BLAS-7 and MTLRS-1 laser systems. These points represent 
the mean value of the full rate data obtained during a 2-minute 
period, and the associated one sigma error bars are indicated on 
each of the 2-minute means. Points which occupy the same 2-minute 
period are differenced (M0BLAS-7 minus MTLRS-1) and plotted as a 
function of the elapsed time. The 2-minute overlap geometry, in 
this example, is designated as 8, 3*, 7. This implies that eight 
consecutive 2-minute differences were obtained, followed by three 
consecutive 2-minute gaps, ending with seven consecutive 2-minute 
differences. Hence a total of 15 minutes of overlap data was taken 
during the pass. The weighted mean difference of these points was 
-1.82 ± 0.11 cm This implies that the MTLRS-1 was measuring short 
with respect to the M0BLAS-7 laser during the pass. A total of 13 
LAGEOS passes was analyzed in a similar manner. 


SIGNIFICANCE 

Results for the 13 LAGEOS collocation passes analyzed are given in 
Table 1. In 12 of the passes the MTLRS-1 laser was found to measure 
short relative to the M0BLAS-7 laser. The weighted mean of the 
entire data set implies MTLRS-1 is short by 1.94 cm. In addition 
to this sytematic bias MTLRS-1 data exhibited systematic jumps of 
as much as 2 centimeters in its range residuals between adjacent 2- 
minute mean points. Analysis for the possible causes of these 
systematics is continuing. 


FUTURE EMPHASIS 

The accuracy of the laser systems being developed continues to 
improve. New technological developments, such as the Micro Channel 
Plate (MCP) will be added to the M0BLAS-7 laser system to improve 
its precision. Preliminary results indicate the rms will improve 
from 2 1/2 centimeters to 1 centimeter. Several collocations 
between laser systems are planned for next year. 
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Table 1. MTLRS-1 Minus MOBLAS-7 Mean Difference 




TWO MINUTE 

WEIGHTED MEAN 

DATE IN 

TIME, 

OVERLAP 

DIFFERENCE, 

1985 

GMT 

GEOMETRY* 

CM 

5/10 

0934 

7 

-1.97+0.21 

5/15 

0612 

12 

1.46±0.14 

5/20 

0628 

12 

-2.00±0.12 

5/21 

0838 

3,1*,6,1*,1 

-0.33±0.15 

5/30 

0704 

6,2*,6 

-1.15±0.10 

6/03 

0848 

11 

-1.52±0.14 

6/14 

0738 

2,1*, 10 

— 2.85 ±0. 1 1 

6/17 

0704 

2,1*, 16 

— 3.80±0. 1 1 

6/26 

0520 

1,1*, 6, 6*, 3 

— 3.77 ±0. 18 

6/27 

0400 

8,3*,7 

— 1 . 82 ±0. 1 1 

6/27 

0732 

18 

— 3.25 ±0. 1 1 

7/04 

0838 

5 

— 1 . 12 ±0.30 

7/14 

0530 

3,1*,3,2*,5 

— 0.29±0.21 


WEIGHTED MEAN = -1.94+0.04 


♦NO OVERLAP 
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Figure 1. LAGEOS Collocation Data 
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GEODYN SYSTEMS DEVELOPMENT 


Barbara H. Putney 
Theodore L. Felsentreger 


OBJECTIVES 

The purpose of the GEODYN Orbit Determination and Parameter 
Estimation, the SOLVE, and the ERODYN computer software is to 
recover and analyze the recovery of geodetic and geophysical para- 
meters from satellite and other data in a state-of-the-art manner 
and in the most operationally efficient manner. 


BACKGROUND 

In 1971 the NONAME and GEOSTAR programs were combined to create the 
GEODYN program. The SOLVE program was created at the same time. A 
few years later the ERODYN, error analysis program was written. 
The philosophy of the development of the software system has been 
maintenance of computer-efficient, well -structured software, with 
appropriate orbit, earth and numerical models, using precise 
satellite measurement modeling and efficient numerical models, and 
performing careful benchmark procedures. This care has paid off in 
the production of several GEM's (Goddard Earth Models), precision 
station locations, improved tidal, GM, polar motion and earth 
rotation values, consistent baselines, and encouraging GRM 
simulations. Careful usage analysis, and modeling using laser, 
Doppler, altimeter and other satellite data from LAGEOS, SEASAT, 
STARLETTE, GEOS and BE-C satellites as well as many others has made 
these accomplishments possible. 


RECENT ACCOMPLISHMENTS AND SIGNIFICANCE 
A. Current Software Usage and Development 

Continued solutions for gravity field, tides, pole positions, 
earth rotation, GM, and baselines, and collocation studies has been 
made as part of the Crustal Dynamics Project. This has been a 
transition year for the research teams from the GEODYN I to the 
GEODYN II program. 


PRECEDING page blank not rlmed 
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Using the GSFC TOPEX Gravity Model effort as the goal, the 
GEODYN II program received its final shakedown, and differences 
between GEODYN I and GEODYN II were reconciled. Careful comparisons 
between preprocessi ng corrections for optical, laser and Doppler 
satellite data were made. Algorithms such as the expanded tidal 
model and old TOLS option (using only part of the geopotential in 
the computation of the variational equations) were modified and 
i ncorporated to shorten the required computer time for the GEODYN 
II program. Detailed timing comparisons were made between GEODYN I 
and GEODYN II on the Cyber 205 and on the IBM 3081 and Amdahl V7 
computers. All-in-all, the software performed to expected timing 
goals and for normal equation generation far better than antici- 
pated, but showed there is still room for further timing 
optimization. This is especially true in the IBM version of GEODYN 
II. 

The SOLVE program, which was vectorized a few years ago, 
continued to perform flexibly to produce the required solutions 
effectively. The new GEODYN II E-matrix format was incorporated 
into the software. Optimizing the I/O (input/output) in the system 
is the logical next step for SOLVE, since on large inversions the 
I/O takes twice as much time as the CPU time. A technique for 
doing concurrent I/O has been employed with some success. Master/ 
slave station constraints and a station network constraints are now 
SOLVE options. Some automatic documentation software have been 
developed for "flowcharting" the program and cross referencing 
the subroutines and common blocks. A continued effort to make the 
program more user-friendly has motivated a change to the program 
input now allowing ranges of parameters to be suppressed and 
adjusted, and in the UCL (Job Control Language) section by having 
the software determine the required size of files and allocating 
the space from within the program. 

The ERODYN program, which was recently converted to the Cyber 
205, has been modified to use the GEODYN II E-matrix format. The 
matrix multiplication and inversion algorithms have been vectorized 
and the I/O has been modified for the Cyber. Orbit error propagation 
has been implemented and tested. The software has been updated to 
include the acceleration parameters now available in GEODYN. 

In addition, there were several specialized capabilities 
requested by users which were incorporated into the ERODYN program. 
For example, the capability of computing the rms effects of unadjus- 
ted parameters and the total rms effects over a whole arc were 
implemented. 
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B. GEODYN II Development (Vector GEODYN) 

GEODYN II has been designed to run in two parts, GEODYN IIS 
(starting) and GEODYN I IE (ending). This was done so that the 
initialization stage of the program, which is sequential in nature, 
would not waste precious Cyber 205 computer time. Currently, GEODYN 
IIS runs only on an IBM compatible front end to the Cyber 205. The 
output of IIS is two files, an interface file and an observation 
file. The GEODYN IIE program can run on either the IBM compatible 
computer or the Cyber 205 Vector Processor. 

The new software has been designed to process a pass of satel- 
lite observations at a time instead of the one observation at a 
time. This design has produced the vectors that are needed for 
efficient utilization of the Cyber 205. This required orbit infor- 
mation be stored for the entire pass of data, which implied a 
rewrite of the software. 

Currently, the program can process laser, optical, and Doppler 
data and all related observation models, which includes biases and 
center-of-mass offsets and required preprocessing corrections. 
Still missing are the ocean-loading, and albedo models. Plate 
motion and GEOPOL (movement of the earth's crust through polar 
motion in the gravity computation) have been incorporated. The DE 
200 JPL Planetary Ephemeris using the Wahr nutation series is now 
usable by the software. 

C. GEODYN II Timings and Their Significance 

Figure 1 shows several satellites and the ratio between GEODYN 
I and GEODYN II Cyber 205 computer charges which is essentially CPU 
time. The runtime ratios vary from 2.5 to 4.2. This is a good 
beginning but it is anticipated the ratios could be made signifi- 
cantly higher with efficiency analysis leading to further 

optimization. The expanded tidal model, not represented in this 
figure, has been optimized already showing significant improvement. 

Figure 2 presents STARLETTE E-MATRIX TIMINGS for various 
numbers of adjusted parameters and stepsizes. The effect of misuse 
of computer memory is demonstrated with the first two entries in 
the figure. By simply adding 12 more large pages (LP is 65000+ 
words) of memory the STU (System Time Unit-like seconds) drops from 
437 to 95. These runs were made with a 2 million words of memory 
which is in the process of being upgraded to 4 million words. The 
Cyber is amazingly fast as is shown by the largest of these runs 
only using 382 STU's. This will allow us to pursue large modeling 
efforts that were impossible in the past. 
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Figure 3 graphically demonstrates the relationship between the 
number of adjusted parameters, STU's required, size of the 
geopotential, and the effect of tide adjustment on two satellite 
arcs, GEOS-1 and STARLETTE. The size of the STU's makes it reason- 
able to do the largest solutions for these satellites. The flat- 
ness of the curve is due to vectorization and memory size. 


FUTURE EMPHASIS 

The final shakedown of GEODYN II has occurred this year for laser, 
optical and Doppler data. At this point the altimeter model and 
satellite-to-satel 1 ite model are partially implemented. Both will 
be completed. Generation of normal points will be added directly 
to the GEODYN program. It is also hoped that some data simulation 
capability will at least be initiated. Adding models, such as 
ocean loading, that bring the GEODYN II program up to the capabili- 
ties of the GEOYDN I program is another important goal. Continued 
emphasis on program efficiency and optimization will mold this 
software into a more useful tool. 

The SOLVE program will continue to evolve to maintain capability 
with the GEODYN program. Further optimization in the arc elimina- 
tion area has become a necessity with the large number of arc 
parameters now being used (combining several arcs maintaining the 
arc parameters). It is still a goal to create a master source file 
containing both IBM and Cyber programs. This is a time consuming 
task and is making some progress, but usually gets interrupted for 
more pressing program modifications. 

The ERODYN program will continue to be evolved towards vectorization 
and compatibility with the new GEODYN. Making the input and output 
easier to understand for the user, is another goal. 

Effort will continue to make the Crustal Dynamics analyses, MARS 
simulation, and TOPEX gravity model determination as successful as 
possible. Software will be modified, optimized and created as is 
required to accomplish these objectives. 
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Timing Comparison: GEODYN I vs. GEODYN II TVM 09/10/85 
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Figure 1. Timing comparison: GEODYN I vs. GEODYN II 


STARLETTE E-MATRIX TIMINGS 
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Figure 2. Starlette E-Matrix timings 



THE APPLICATION OF COMPUTER IMAGING TECHNIQUES TO 
OCEAN TOPOGRAPHY AND GEOPHYSICAL DATA 


Carole Nethery 


OBJECTIVE 

Understanding of scientific data may be improved by converting it 
into an image that may be studied visually. The scientist can 
identify features, patterns and common elements that would otherwise 
go unnoticed. The analysis may be done interactively or by viewing 
photographic copies of the image. Photographs of images provide a 
means for the scientist to communicate findings to others in his/her 
professional community. The purpose of this particular set of 
image analysis programs is to extend the knowledge of the geophysi- 
cal properties of the oceans. This aim assumes the use of the 
advanced techniques available at the Goddard Space Flight Center. 
Mindful of the merits of software sharing and the importance of 
developing interdisciplinary software, it is also assumed that 
programs will be designed with generality so that they may be used 
for a variety of applications. 


BACKGROUND 

In 1982, several programs were written to convert SEASAT altimeter 
data describing the ocean surface to images. The software design 
objectives were: efficient computer use, ease of modification and 

program generality. After several changes in hardware and user 
personnel, isolation of hardware-dependent code and facility of 
user input were given equal consideration. The latter was carried 
out in part by modifying programs to operate under TAE, 
(Transportable Applications Executive). Interfaces to the TAE 
catalog manager were written to simplify user-file interaction. 

The programs written handled tape-to-disk transfer, image mapping 
that allows screening of areas and data ranges, image display, 
pseudo-coloring of images, calculation of image statistics to aid 
coloring, and simulation of the illumination of an image by a 
rotating light source. Procedures were developed for photographing 
images for publicaiton. The resulting images provided information 
for publicaitons on mean sea surface and mesoscale ocean surface 
variability studies. 
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RECENT ACCOMPLISHMENTS 


Mean sea surface data, bathymetry data and gravity field models 
computed from satellite altimetry have been processed to form 
photographic images. An example of illuminated topography based on 
bathymetry data is shown in Figure 1. This image was based upon 
the 5' x 5‘ grid of data provided by Washington University in St. 
Louis. The image will be prepared in poster form as a companion to 
the mean sea surface maps. 

Software has attained production status on the LTP Vax. The ver- 
sions of the programs installed have improved running times, are 
easier to use and have produced more consistent photographic genera- 
tion than in the past. The image terminal interface isolates 

hardware dependent code. A reasonably efficient procedure for 
image restoration has been established. Programs have been written 
to merge files and subset large global files. The tape-to-disk 
program to read data runs more efficiently. Efforts have been made 
to make the software easier to use by other applications. 

Color enhancement has been improved, particularly for photographic 
output. The photographic software colors data for the Optronics 
hardware. Since the eye can see literally millions of colors, 
color enhancement is complicated and the possible results powerful. 
The CIE (International Commission of Illumination) System, will be 
tried as a base. This system describes colors mathemtical ly and is 
founded on measurements of color-matching abilities of the average 
eye. This more logical system might make it easier to handle 
discrepancies between calculated, screen, and photographic color. 

New software is being tested to manage files. Problems with the 
catalog manager of TAE and unnecessary running time lost has promp- 
ted this effort. The management will tailor the file handling to 
our needs and allow us to bypass those items that we do not use. 
Groupings of files will be supported and informative listing and 
selections provided. It is planned that the file selection proce- 
dure be available to the IBM PC's. Files foreign to the ocean 
system, yet appropriate for processing, will be accepted with the 
user suppling the necessary parameters. 
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SIGNIFICANCE AND FUTURE EMPHASIS 


The programs have provided oceanographers and geophysicists with 
new information about the oceans and broadened their understanding 
of the data. Images have appeared in publications and have been 
presented at scientific meetings and conferences. The image of the 
global mean sea surface appears in U.S. News and World Report. The 
image of mesoscale ocean height variability was used for an 
internationally-distributed NASA Poster. The direction of software 
development will continue to follow the needs of the scientists in 
support of their objectives. 

The future emphasis will be on image and graphic processing of 
multiple data sets. Design will lean toward the use of personal 
computers for software development, with efficient links to the VAX 
for image display, and MPP (Massive Parallel Processor). The MPP 
which can handle many calculations at a very high rate, and whose 
design adapts itself to image processing is almost ready for 
application. We will be using the MPP, as it becomes available for 
most image preparation work, remapping and illumination. The new 
Merlin three-dimensional graphics terminal will be connected to the 
MPP as well. The Optronics photographic procedure will probably be 
used initially for the graphics output. 
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Figure 1. Sea surface topography based on bathymetry. 


GEODYNAMICS LASER RANGING SYSTEM FOR AN EARTH OBSERVING 


SYSTEM SATELLITE 
Steven C. Cohen 


BACKGROUND 

NASA is currently studying the concept of an Earth Observing System 
(Eos) (NASA, 1984) which involves multisensor satellites making a 
variety of Earth observations from platforms in near-polar orbits. 
One of the planned instruments for inclusion on an Eos satellite is 
the Geodynamics Laser Ranging System (GLRS) (NASA, 1986) GLRS 
consists of a high precision, low energy, low mass NdrYAG laser, a 
precise pointing system, a fast response photodetector and a 
variety of associated optoelectronic components. GLRS would be 
used to make very accurate geodetic measurements by ranging to cube 
corner retroreflectors on the Earth. It could also be used to make 
land, ice, and ocean altimetric measurements by ranging at nadar. 
The laser ranging measurements would be made to arrays of targets 
with nearest neighbors separated by 10-100 km. The measurements 
would be used to determine target positions and intersite distances 
and to study regional and local scale crustal movements including 
strain accumulation and release in seismic zones, strain propagation 
from plate boundaries into .plate interiors, uplift and subsidence 
associated with both tectonic and nontectonic processes, and 
regional scale intra- and interplate motions. The positions of the 
GLRS ground sites could be tied to those of a fiducial network 
(e.g., the ground-based observatories operated by NASA as part of 
the Crustal Dynamics Project) to provide a global geodetic' network . 
With targets located on ice sheets, the GLRS measurements could be 
used to monitor ice velocities and strain. The GLRS measurements 
would also be used to provide a highly accurate satellite ephemeris 
for Eos to be used by other precision position devices such as a 
microwave altimeter. When operated as a short wavelength altimeter 
GLRS could map the topography of ice sheets and various land forms 
with higher spatial resolution and accuracy than possible with 
other techniques. These measurements would be useful to studies of 
ice dynamics, geomorphology, tectonic plate flexure, and hydrology. 
Simulation studies have begun to assess the sensitivity of the 
deduced geodetic parameters from the GLRS ranging measurements 
(positions, intersite distance, orbital parameters) to noise in the 
observations and uncertainties in model parameters. 
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RECENT ACCOMPLISHMENTS AND SIGNIFICANCE 


Error analysis and mission simulation studies have begun using 
GSFC's GEODYN, ERODYN, and SOLVE software programs. The initial 
studies have focused on determining the precision and accuracy in 
the deduced intersite distances for an array of cube corner targets 
located in California as shown in Figure 1. The targets are located 
every 50 kilometers throughout the state. The parameters used in 
the study are shown in Table I. Most important are the orbital 
parameters including an altitude of 800 km (circular orbit) and an 
inclination of 100 degrees. The laser system has a single shot 
precision of 1 cm and is pulsed at 10 Hz. Typical results for a 
3-day observation period are shown in Table II. Observations are 
made over 11 orbital passes, each pass being 5-10 minutes long. 
The orbits are divided into six arcs: five arcs involve two succes- 

sive passes over the grid and more than one complete revolution of 
the Earth, while one arc involves a single short pass over the 
grid. The predicted uncertainties in intersite distances due to 
laser noise are about 0.2 cm for distances up to 500 km. Predicted 
uncertainties in the vertical are also well less than 1 centimeter. 
The absolute accuracy in the deduced parameters depends on 
the accuracy of the model parameters used in the data reduction. 
The preliminary analysis has focused on the effects of uncertain- 
ties in gravity field, drag, and solar radiation parameters. 
Results show that the gravity field uncertainties are dominant. 
Using covariances that are one-quarter that available from the GEM 
10B gravity model (a conservative estimate for the Eos era) results 
in aliases of less than 1 cm (one sigma) for distances up to 100 km 
and about 2 cm at 500 kilometers. Aliases in the vertical are more 
severe and can exceed 10 centimeters at 500 km. Further calcula- 
tions recently obtained suggest that the alias in the vertical 
measurements (and the cross baseline distance changes) can be 
reduced by tying GLRS measurements to fiducial points defined by 
the Crustal Dynamics network. 


FUTURE EMPHASIS 

A number of effects must be considered to complete the analysis of 
the GLRS system. Variations in the atmospheric refractive index 
due to gradients in temperature and pressure can introduce errors 
into the range measurements if they are not taken into account. 
The effects of uncertainties in tidal parameters, polar motion, and 
other geophysical parameters need to be assessed as do the effects 
of timing and ranging biases. The utility of ties to fiducial points 
will be considered in more detail. The current error analyses will 
be augmented by mission studies involving the processing of 
simulated data. Error analyses and mission simulations for alti- 
metric measurements may also be undertaken. 
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Table I. Simulation Conditions 


ORBIT AND MISSION PARAMETERS: 

ALTITUDE: 800 KM 
INCLINATION: 100 DEGREES 
ECCENTRICITY: 0 

TARGET ARRAY: 157 TARGETS (50 KM SPACING) 

SURVEY INTERVAL: 3 DAYS 

11 DATA ACQUISITION ORBITS (20 DEGREE CUTOFF) 
6 ARCS 

APRIORI SATELLITE INFORMATION 
POSITION: 10 METER (X, Y, Z) 

VELOCITY: 0.01 METER/SEC (XDOT, YDOT, ZDOT) 

APRIORI TARGET LOCATIONS: 1 METER (X, Y, Z) 

SYSTEM PARAMETERS 

RANGE PRECISION: 1 CM 
PULSE RATE: 10 PPS 
ACQUISITION TIME: 10 SECONDS 
DWELL TIME: 0.5 SECONDS 
SLEW TIME: 0.5 SECONDS 

OTHER PARAMETERS 

REFERENCE GRAVITY FIELD: GEM10B 

SOLAR RADIATION PRESSURE AND DRAG INCLUDED 


Table II. Error Analysis — Baseline Noise & Alias 


(cm) 


DISTANCE 

STATION # 

NOISE 

ALIAS 

NORTH (km) 




50 

72 

0.2 

0.7 

100 

65 

0.2 

0.6 

150 

59 

0.2 

0.7 

200 

53 

0.2 

0.9 

250 

47 

0.2 

0.6 

300 

41 

0.2 

0.9 

350 

35 

0.2 

0.9 

400 

28 

0.2 

1.3 

450 

21 

0.2 

1.4 

500 

14 

0.1 

1.6 

(EAST) 



0.4 

50 

80 

0.2 

100 

81 

0.2 

1.0 

150 

82 

0.2 

1.1 

200 

83 

0.2 

1.5 

250 

84 

0.2 

2.0 
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MARK III VLBI TECHNOLOGY: AN EVALUATION OF THE NEW GENERATION 

WATER VAPOR RADIOMETERS 


G.L. Lundqvist 
M.W. Hayes 
T.A. Clark 


OBJECTIVE 

The objective of this task is to evaluate the new generation of 
water vapor radiometers (WVRs), by data acquisition and analysis, 
as being capable of measuring the excess microwave path delay due 
to atmospheric water vapor. 


BACKGROUND 

The NASA Crustal Dynamics Project has a total of five water vapor 
radiometers (WVRs) in operation to support VLBI measurements. 
Four of these are an older JPL design by G, Resch dating from the 
mid-1970's (the RO series instruments) which have been extensively 
re-worked by the WVR team at NASA/GSFC, Interferometries, Inc., and 
Bendix Field Engineering Corporation. The fifth is the initial 
prototype of a second-generation instrument developed by M. Janssen 
and his co-workers at JPL (called the JO series). 

The WVRs measure the equivalent sky brightness temperature at 20.7 
and 31.4 GHz. The lower frequency is close to the resonant water 
vapor emission line at 22.235 GHz and is therefore sensitivie to 
the integrated amount of water vapor along the line of sight. The 
second frequency is used to measure the amount of liquid water 
along the same line of sight. By linearly combining the brightness 
temperatures at these two frequencies, when weighted by properly 
chosen retrieval coefficients, the microwave excess path delay 
along the line of sight can be determined. 

Typically, the total delay through the troposphere is on the order 
of 2 meters. Of this, the largest portion is due to the refractive 
effects of nitrogen and oxygen; this "dry" part can be calibrated 
by using an accurate barometer to "weigh" the total number of 
molecules overhead. While the contributions of water vapor are 
much smal ler--typical ly only about 10% of the total— the variability 
(both with time and in direction) is quite severe. Since the CDP's 
precise geodetic measurement programs aim at 1 cm accuracy, the 
WVRs will play an important part in reaching that goal. 
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Traditionally the wet part of the excess delay has been determined 
by using atmospheric models and surface meteorological data. Since 
the correlation is often very poor between surface and higher layer 
humidity the surface model tends to give erroneous values for the 
wet path delay. 


RECENT DEVELOPMENTS 

The CDP's Mojave Base Station has had one of the upgraded RO-series 
WVRs, R07, in operation since December 1984. With this station 
being located in the California desert, the content of water vapor 
can be quite low and at times shows little variability from day-to 
-day. The WVR at Mojave has shown very good long term stability 
with calibration accuracy typically in the 0.3% range (this 
reflects an order-of-magnitude improvement over the performance 
before the upgrade). 

As of March 1986 the VLBI facility at the STDN station on the island 
of Kauai, Hawaii, has the R05 WVR on line. This site has a very 
interesting local topography. To the east is a mountainous rain 
forest with the highest annual rainfall on earth; in the opposite 
direction there exist semi-desert conditions typical of a tropical 
island. This situation should produce a strong azimuthal dependence 
in the data. 

This is in fact what has been found in the data acquired at this 
site. Figure 1 shows an example of how quickly the excess path 
delay can vary although sky conditions and surface humidity showed 
very little indication of variations during the same period. The 
upper part of the plot shows the sky brightness temperature (in K) 
in channel 1 as a function of time and the middle plot shows the 
corresponding data for channel 2. The lower part is the excess 
path delay (in cm) as derived from the two channels above. It can 
be seen that the path delay drops from about 10 cm down to 3-4 cm 
in only 3 hours and then increases again. Figure 2 shows an example 
of the azimuthal dependence. The WVR has measured the path delay 
at 30 degree intervals around the horizon at a constant elevation. 
The variations in path delay as a function of azimuth can be clearly 
seen. The variations have a peak value of 4-5 cm. This variation 
is extreme among the VLBI stations examined so far indicating that 
the WVR is a necessary part of the Mark III system at that site. 

During a recent side-by-side comparison of the R07 and the J01 WVRs 
at the Mojave Base Station some interesting results were found. 
The two WVRs followed the same measurement program on the sky and 
the output of the two insruments was compared. Figure 3 shows the 
results of this comparison. The two WVRs measured the sky bright- 
ness temperature to within 0.2 K which is a remarkably good result. 
The two i ntsruments use very different techniques for the measure- 
ments and this indicates that the absolute accuracy of the sky 
brightness temperatures is on the order of 0.5 K which is close to 
the goal of the instruments. 
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SIGNIFICANCE 


The significance of these results is that the instruments seem to be 
capable of measuring excess wet path delay to an accuracy of 0.5 - 
1.0 cm. The VLSI data quality has been improved by as much as a 
factor of 2 on single experiments using path delay corrections 
from the WVRs. 


FUTURE EMPHASIS 

By the summer of 1986 the CDP should have WVR's at both Westford 
and Fort Davis which should provide some very powerful data to 
apply to the VLBI experiment data sets in order to verify the 
usefulness of the WVR instruments. 
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HIGH RESOLUTION GRAVITY FIELD MAPPING 
WITH A SPACEBORNE GRAVITY GRADIOMETER 


Werner D. Kahn 


OBJECTIVE 

The gravity field of the earth reflects variations in the earth's 
composition. Knowledge of the earth's gravity field provides 
information on the density structure of the earth's interior and 
outer layers from which we can increase our understanding of the 
formation of geological features, such as mountains and basins, and 
the mechanisms and scale of processes such as mantle convection. 
In the area of ocean science, the gravity field is used to construct 
the reference surface (geoid) for determining the dynamics of the 
global sea surface. An ultrasensitive gravity gradiometer seems to 
be a suitable sensor for mapping the fine structure of the earth's 
gravity field. Studies have thus been undertaken to assess that 
instrument's capabilities. 


BACKGROUND 

Covariance error analysis techniques have been used to assess the 
capabilities of a spaceborne gravity gradiometer to map the fine 
structure of the earth's gravity field. Since the gradiometer 
measures the second derivatives of the geopotential field, the 
gradiometer will be sensitive to the density variations of the 
earth's outer crust which is the source of the high frequency 
components of the earth's gravity field. A gravity gradiometer 

mission would provide in situ measurements on a global scale of the 
earth's gravity field fine structure. For the covariance error 
analyses, the effects of orbit altitude, orbit accuracy, and gradio- 
meter precision have been examined. 


RECENT ACCOMPLISHMENTS 

A computer program has been developed and tested to perform error 
analysis studies to assess the capabilities of a spaceborne gravity 
gradiometer. The program was designed to recover gravity or geoid 
anomalies for the full gravity gradient tensor, the three axes, 
and the single axis gradiometer configuration. Error sources for 
the gradiometer studies were instrument precision, orbit position 
accuracy, satellite attitude and attitude rates. The program 
provides estimates of high frequency components of the gravity 
field over a local region. Both mean square gravity and geoid 
anomalies as the solution parameters were estimated by this program. 


PRECEDING PAGE PLANK NOT FILMED 

209 


The mathematical model for the disturbance potential was derived 
from Poisson's integral (Heiskanen and Moritz). The components of 
the gravity tensor were developed by forming the second order 
parti als from the geopotential and a sequential estimator was used 
to provide estimates of the gravity and geoid anomaly uncertainties 
based on the modeled error sources. A detailed report describing 
the mathematical analysis is in preparation. Results of the error 
analyses to be discussed below assumed the 3-axes gradiometer 
configuration which in mathematical terms are the diagonal elements 
of the gravity tensor. They satisfy Laplace's equation (i.e., 
V 2 T=0) and only two of the three components are independent. The 
3-axes gradiometer is the configuration now under development by 
Professor H.J. Paik at the University of Maryland. 

Figure 1 shows the evolution of gravity and geoid anomaly uncertain- 
ties as a function of resolution. The results depicted are for a 
orbiting gradiometer in a 160 km and 200 km polar/circular orbit 
and a mission duration of 180 days. The gradiometer precision of 
10“ 3 EOTVOS Units (E)* with, a measurement sampling rate of 1 mea- 
surement per 4 s , and a priori knowledge of the gravity field of 9 
milligals/l 0 blk and 50cm/l° blk respectively were also assumed . On 
the average there are approximately 60 measurements in each l°xl° 
blk which corresponds to 16 orbital passes over the block in a 180 
day period. Under these assumptions, it can be seen that for the 
satellite in a 160 km altitude orbit, 1° x 1° gravity and geoid 
anomalies can be estimated with an accuracy of 0.23 mgals and 1.2 
cm respectively. At an orbital altitude of 200 km the above results 
are 0.45 mgls and 3.5 cm respectively. It can also be seen that 
the degradation in accuracy is more sensitive to geoid anomalies 
than to gravity anomalies when increasing the orbit altitude from 
160 km to 200 km. The accuracy is thus degraded by a factor of 2 
for gravity anomalies and a factor of 3 for Geoid anomalies. The 
results of Figure 1 show the 3-axes gradiometer to be well suited 
for mapping the high frequency components of the geopotential. 
Even at a 200 km orbital altitude, the global gravity field can be 
very accurately mapped with an orbiting gravity gradiometer having 
10~ 3 E system precision. 

In Figure 2 the effects of gravity gradiometer system precision is 
examined for orbit altitudes of 160 km and 200 km, respectively. 
The results are presented in terms of l°xl° gravity anomaly accuracy 
only. It can be seen that the same levels in accuracy due to 
increase in orbit altitude prevail as those presented in Figure 1. 
It is interesting to note however, no significant increase in the 
accuracy of estimated gravity anomalies results from increased 
gravity gradiometer precision. The curves are relatively flat when 
the precision in the gradiometer increases from 5xlO" 3 E to 10“^E. 
Large changes do however occur when the precision is increased in 
the range 10 _1 E to 10 _2 E. These results indicate that a spaceborne 
gravity gradiometer at an orbit altitude of 160 km to 200 km can 


* IE = 10- 9 /s 2 
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with a precision of 10"3 e to 6xlO“^E, map the high frequency compo- 
nents of the earth's gravity field to an accuracy better than 1 mgal 
in terms of gravity anomalies and better than 4 cm in terms of 
geoid anomalies with a horizontal resolution of 110 km. 


FUTURE EMPHASIS 

The computer software will be augmented to include the satellite to 
satellite tracking (SST) measurement data in combination with 
gravity gradiometer data. Studies will then be performed to eval- 
uate the level of improved accuracy in gravity and geoid anomaly 
estimation resulting from the combination of these data types. 
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Figure 2. Estimation of 1 0 X 1 0 Gravity Anomalies vs. Gradiometer Precision 
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THE POWER CURVE RECOVERY OF DIURNAL TEMPERATURE EXTREMA 


Jean E. Welker 


OBJECTIVES 

The objectives of this work are: 1) To describe and model the power 
curve recovery of diurnal temperature extrema, DTE. 2) To discrimi- 
nate between the different types of precipitation cooling events 
which cause drops in the ( diurnal temperature maxima, DTMAX, in 
soils, air, and shallow pan evaporation water. 


BACKGROUND 

From a one-dimensional view of temperature alone, variations at the 
earth's surface manifest themselves in two cyclic patterns of 
diurnal and annual periods, due principally to the effects of 

diurnal and seasonal changes in solar heating as well as the gains 
and losses of available moisture. A related but very different 
third cycle has been identified, with a variable mesoscale period 
of from 3 to 14 days. This is the period of recovery or return to 

the approximate level of original values of a series of 3 to 14 

monotonical ly increasing DTE variables, at 10 cm depth below bare 
soil after a precipitation-cooling event. An idealized schematic 
of such an individual cycle is shown in Figure 1, where the 

recovery leg of 3 to 14 days is represented by 2 , which begins at 
the minimum value of a. These series of values are obtained from 
raw data records; a minimum of three consecutive monotonical ly 
increasing values is a requirement for these restricted series of 
values. This third cycle is dependent on the variables of diurnal 
temperature extrema, DTE, extreme temperature values recorded each 
day at a test site; the diurnal temperature cycle, on the other 
hand, is related to the very different variable of temperature 
readings taken at fixed intervals throughout the day. Although the 
daily temperature high occurs many times in the afternoon, and the 
daily temperature low occurs many times shortly after midnight, the 
DTE variables are not constrained absolutely by a specific time of 
day, and only are affected by the diurnal cycle in an integrated 
way. These variables do undergo an annual period by rising in 
value in the summer and declining in value in the winter. The DTE 
variables in their dual forms, diurnal temperature maxima, DTMAX, 
and diurnal temperature minima, DTMIN, each have the peculiar 
characteristic of recovering from a precipitation-cooling event in 
a power curve, when the raw data series are a sequence of monotoni- 
cally increasing values. The DTE recovery clearly depends on solar 
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heating of the soil with an increased moisture content from pre- 
cipitation combined with evaporation cooling at soil temperature 
lowered by precipitation cooling, and many other variables, but 
is quite regular and universal for vastly different geographical 
locations, and soil types and structures. The DTE variables 
gradually rise in value through the summer months from the effects 
of the alternate heating and cooling of the diurnal period, as seen 
from the raw data values for the DTMAX variable in Figure 2, for 
the three soil depths of 5, 10 and 20 cm. This slowly changing 
growth of DTMAX in a quasi-equilibrium condition is frequently 
interrupted by a precipitation-cooling event. When a cooling or 
drop in the value of DTMAX occurs, the initial increments of change 
are large as the variable recovers and approaches the equilibrium 
value it had prior to the cooling event. The characteristics of 
an ascending power curve are also initial large increments of 
change which become progressively smaller. This recovery occurs 
slowly over the mesoscale period of 3 to 14 days, so long as 
another precipitation cooling event does not occur in the meantime, 
which would restart the recovery cycle after a new minimum DTE 
value had been achieved. The regularity of the power curve 

recovery has allowed an accurate power curve fit to the raw DTE 
data, i.e., each recovery leg, 2, shown in Figure 1, has been 
fitted with a power curve, for a value of the minimum value, a, and 
a value for the power exponent, b. From these data fits, a predic- 
tive model approach which generalizes the DTE recovery for all 
precipitation-cooling events has been developed. The generality of 
the DTMAX power curve recovery for soils was extended by the com- 
parison of its properties with respect to water in evaporation pans 
and air at the same test sites. 

Although the research on the manner of recovery of the DTMAX variable 
has progressed, a more complete interpretation of the process 
requires a greater understanding of the initial precipitation 
cooling process. A precipitation event which cools the soil can 
affect many variables relating to soil temperature, and the 
characteristics of the precipitation event are, in themselves, 
variables. Much of a heavy and intense rainfall over a short dura- 
tion can run off and not be absorbed in the soil; a slower, longer 
duration rainfall can have more of its total moisture absorbed in 
the soil, under the same conditions. The time of the precipitation 
event is also important, for summer afternoon precipitation on hot 
soil will differ in effect from the same amount of precipitation 
on cooler soil in the night. In all cases, the precipitation event 
occurs over some duration of time relative to the time at which the 
value of DTMAX in the soil occurs; under ideal conditions, DTMAX 
occurs sometime in the afternoon. 
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RECENT ACCOMPLISHMENTS 

Recent accomplishments are three-fold: the characteristics of the 
power curve recovery of the DTE variables are now being considered 
for soil at different depths, air, and water in evaporation pans. 
The fact that the power curve recovery occurs in all these sub- 
stances provides insights as to the energy balance conditions which 
cause this recovery characteristic. A number of data sets have 
been entered, analyzed and plotted and have supported initial 
conclusions. A second approach to the problem is to investigate 
energy balance experiments which have been conducted for a number 
of days in order to explain the power curve recovery. One famous 
experiment conducted in Arizona in the early sixties tracked the 
energy balance budget both in soil with changing moisture content 
as well as shallow water and listed recorded value of the DTE 
variables. For this single example, the reconstructed DTMAX series 
followed the power curve recovery of the many other examples tested, 
with the distinct advantage of a detailed description of the energy 
balance components. (Fritschen and van Bavel , 1962). A third 
approach has been to assess more detailed precipitation data in 
order to understand the precipitation cooling process which pro- 
duced the initial drop in the DTE variables. Precipitation cooling 
has been perplexing in the past because of a double-valuedness in 
the amount of precipitation that produces the same drops in the raw 
DTMAX values. To investigate the problem further, hourly precipita- 
tion data have been acquired for one of the station data sets 
previously considered. After much analysis of the data, it appears 
that all precipitation data in the 24 hours prior to the recorded 
drop in the values of DTMAX must be included, and that precipita- 
tion values less than or greater than 2.5 cm can produce double- 
valuedness in the DTMAX values as shown in Figure 3. The implica- 
tions of these results are being investigated. 


SIGNIFICANCE 

The regularity of the power curve recovery of the DTMAX variable in 
soils at different depths, air, and water in evaporation pans 
allows predictive modeling of the recovery as well as affording 
an insight into the relevance and the importance of some of the 
energy balance components and how they affect the DTE variables. 
Because of the frequency of the DTE drops and recoveries combined 
with the duration of the mesoscale recovery period, the DTE vari- 
ables are in the recovery phase a significant portion of the time. 
For the exceptional data series from the Tifton Experimental 
Station in Georgia in 1979, the DTMAX variable was in the recovery 
phase for 69 percent of the March through August period. The power 
curve recovery is so consistent and regular during this phase that 
it can be modeled and predicted with varying degrees of accuracy 
for a number of meteorological stations under a wide variety of 
conditions. 
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FUTURE EMPHASIS 


The energy balance conditions for soils at different depths, shallow 
water, and air will continue to be analyzed and appropriate models 
will be developed. The double-valuedness of the precipitation 
cooling events have only been characterized for very limited data 
sets and will be expanded to include others. Further statistical 
analysis is required to generalize the power curve recovery for a 
larger number of stations and to further investigate the behavior 
of the DTMIN variable. A journal article on a portion of the results 
to date is in preparation. These results had been developed over a 
period of time in consultation with Professor Helmut Landsberg of 
the University of Maryland, whose untimely passing in December 
1985, is mourned by this author and countless others who have been 
helped and befriended by this exceptional person. 
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